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2.  Objectives 

The  objective  of  this  work  is  to  develop  constitutive  models  that  are  capable  of 
accounting  for  the  evolution  of  microstructure  in  finite-deformation  processes  of 
heterogeneous  materials.  The  models  to  be  developed,  through  rigorous 
homogenization  procedures  for  random  microstructures,  will  take  the  form  of 
standard  homogenized  constitutive  models,  which  will  depend  explicitly  on 
appropriate  internal  variables  (serving  to  characterize  the  current  state  of  the 
microstructure),  and  will  be  supplemented  by  differential  evolution  laws  for  these 
internal  variables.  The  particular  case  of  a  porous  metal  will  be  considered  in  some 
detail  and  the  newly  developed  constitutive  models  will  be  implemented  in  a 
general-purpose,  finite-strain  finite  element  program.  Several  problems  in  the  areas 
of  metal  forming  will  then  be  analyzed,  and  implications  for  the  onset  of  shear 
instabilities  will  be  explored.  Finally,  comparisons  with  experimental  results  and 
further  refinements  of  the  models  will  be  carried  out,  for  porous  aluminum  and 
aluminum-matrix  composites. 


3.  Status  of  effort 

This  project  has  been  completed  as  of  12/31/99. 


3.  Accomplishments 

The  developed  constitutive  models  are  the  first  of  their  type  to  be  able  to  account  for 
deformation-induced  anisotropy  in  porous  materials  and  particle-reinforced 
composites.  The  main  findings  of  the  work  are:  (a)  the  evolution  of  the 
microstructure  affects  quite  significantly  the  macroscopic  response  of  porous 
materials;  (b)  there  significant  coupling  between  the  various  microstructural 
variables,  in  particular,  between  the  porosity  and  anisotropy;  (c)  qualitative 
agreement  has  been  found  with  available  experimental  results,  particularly  for  the 
evolution  of  porosity  and  anisotropy  in  forming  processes.  A  constitutive 
subroutine  (UMAT)  for  the  model  has  been  successfully  implemented  in  ABAQUS 
and  tested  for  several  different  t5q3es  of  forming  processes,  including  plane-strain 
and  axisymmetric  extrusion  and  tapered-disk  compaction.  In  the  context  of  future 
work,  we  will  explore  further  the  implications  of  microstructure  evolution  on  the 
overall  stability  of  deformation  processes  and  consider  applications  to  problems 
where  microstructure  evolution  is  significant,  including  ductile  fracture,  high-strain 
rate  (explosive  loading)  applications  and  impact  loading. 
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project.  In  addition,  one  former  graduate  student,  M.  Kailasam,  now  employed  by 
HKS  (the  makers  of  ABAQUS),  assisted  with  the  computations  on  his  free  time.  In 
addition.  Professor  Nick  Aravas,  now  at  the  University  of  Thessaly  (Greece) 
collaborated  with  us  on  the  numerical  implementation  of  the  constitutive  model 
and  in  the  writing  of  the  papers. 
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Transitions 


Re-initiated  contact  with  Dr.  Richard  Becker,  formerly  at  ALCOA,  and  now  at 
Lawrence  Livermore.  Dr.  Becker  has  proposed  the  implementation  of  our 
anisotropic  constitutive  model  for  porous  metals  into  the  computer  programs  at 
LLNL,  as  part  of  the  ASCI  initiative.  If  the  project  is  approved,  the  constitutive 
model  will  be  extended  for  the  materials  of  interest  at  LLNL  and  implemented 
numerically. 
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in  Solid  Mechanics  and  Materials  Science,  University  of  Cambridge,  September- 
December  1999. 

Invited  to  participate  and  lecture  in  the  Workshop  on  "Homogenization  and 
Effective  Media  Theories."  Mathematical  Sciences  Research  Institute,  Berkeley, 
March  6-17,  2000. 

Invited  to  give  a  "Sectional  Lecture"  on  "Nonlinear  Composites"  at  the  upcoming 
20*  International  Congress  of  Theoretical  and  Applied  Mechanics,  August  2000. 


9.  Appendix:  The  finite  element  program 


o  o 


l^-f 

-  UMAT  routine 

subroutine  uinat (stress, statev, ddsdde, sse, spd, scd, 

&  rpl,ddsddt,drplde,drpldt, 

&  stran,dstran, time, dtime, temp, dtemp,predef , dpred, cmname, 

&  ndi,nshr, ntens,nstatv,props, nprops, coords , drot, pnewdt, 

&  celent,  dfgrdO,dfgrdl,noel,npt, layer, kspt,kstep,kinc) 

include  ' ABA_PARAM . INC ' 
character* 8  cmname 

real*8  stress (ntens) , statev (nstatv) , ddsdde (ntens,n tens) 
real *8  ddsddt (ntens) ,drplde (ntens) , stran( ntens) ,dstran (ntens) 
real *8  time (2) ,predef (1) , dpred (1) , props (nprops) , coords (3) 
real*8  drot(3,3) , dfgrdO (3 , 3) , dfgrdl (3 , 3) 
real* 8  sse, spd, scd, rpl , drpldt, dtime, temp, dtemp 
real * 8  pnewdt , celent , ct , s t , cp , sp , cs , ss 
integer  ntens , ndi , nshr, nstatv, nprops , noel , npt , layer 
integer  kspt, kstep, kinc, i, j , inpt 

real *8  mule, kle,mu2e, ek2e 

real* 8  s(3,3) , strain (3, 3) , ddsdde_mine ( 6 , 6 ) 
real* 8  f , wil , wi2 , dlam, atheta, aphi , apsi 
real*8  sc (6) ,strainc(6) ,dummy3(3,3) 
real*  8  dummyl(3,3)  ,R(3,3),RT(3,3)  ,  duinmy2  ( 3 , 3 ) 
real * 8  mul , kl , mu2 , k2 , sigyl , x2 , xacc , s igyO , xl , hrd 
real *8  eps_plastic, kll, rtsafe, rot (3,3), rott (3,3) 
real* 8  rtnewt, eigen (3, 3) ,ei2 (3) ,ei3 (3) ,eil (3) 
real *8 . a, b, c, xbl, xb2 
real *8  func 
real *8  eigval(3) 
real*8  rota (3, 3) , rotj (3, 3) 
real*8  tottheta, totphi, totpsi 
integer  nb, num, n, np, nrot 
logical  path, yield, check, neg, debug 
logical  evolf , evolwi , evolti 
logical  zb, cal 

common  /controls/evolf , evolwi , evolti 
common  /mkmodulus  /mul ,  kl ,  mu2 ,  k2  ,  s  igy 0 ,  kl  1 
common  /mkelas/mule,  kle,mu2e,  ek2e 

common  /mkdatal/f, wil, wi2 , sc, strainc, eps_plas tic, sigyl, hrd 
common  /mkorient/cp, sp, ct, st, cs, ss, rot, rott, atheta, aphi, apsi 
common  /mkdata2 /yield, checks neg 
common  /mkrot/R, RT, rot j 
common  /mdebug/ debug 
common  /call/cal 
common  /zbl/zb 

external  func, rtsafe, funcd, rtnewt, zbrak 
path=: .  true. 

open (unit=15 , f ile= ' / scratch/mahesh/ data ' , 

&status= 'unknown' ) 

open (unit=16 , f ile= ' /scratch/mahesh/datal ' , 

&status=' unknown' ) 
c 

Q -  Various  parameters,  like  those  used  to  turn  the  evolution  of 

c  some  of  the  micro.  varicUDles  on/off  etc.  are  entered  below, 

c  if  ( ( (noel.eq.150) .or.  (noel .eq. 600) ) .and. (npt . eq. 1) )  then 

write(15, *) 'kinc' ,kinc, 'npt' , npt, 'ntens ', ntens 
writedS,  *)  'stress' 
write (15, *) stress 
c  cal=.true. 

call  flushdS) 
c  else 

c  cal=. false, 

c  endif 

evolf = . true . 

evolwi =. true. 

evolti=. true. 

debug= . false. 

if (debug. eq. . true . )  then 

write (15, *) ' entry' 

call  flushd5) 

endif 

inpt=npt 

yield= . false. 

neg= . false . 

mul=props (1) 

kl=props (2) 

mu2=props (3) 

k2 nprops (4) 

sigyO=props (5) 

mule=props (6) 

kle=props (7) 

mu2e=props (8) 

ek2e=props (9) 


kll=kl*1.0 

if (debug . eq . . true . )  then 
■  write(15, *) 'elas ' ,mule, kle,mu2e, ek2e 
call  flush (15) 
endif 

do  1  i=l,3 
do  2  j=l,3 

if  (dabs(dfgrdO(i, j) ) .lt.l.Od-14)  then 
dfgrdO (i, j ) =0 . 0 
endif 

if  (dabs (dfgrdl (i , j ) ) . It . 1 . Od-14)  then 
dfgrdKi,  j)=0.0 
endif 

continue 
continue 
do  3  i=l,6 

if  (dabs(stress(i) ) .lt.l.Od-08)  then 
stress (i) =0.0 
endif 
continue 

Convert  stress  to  column  form  in  MY  notation 

The  input  comes  in  the  form:  {sll / s22 # s33 , sl2 / sl3 / s23 } 

while  I  have  always  used 

{sll,s22,s33,sqrt(2.0) *sl2 , sqrt (2 . 0) *s23 , sqrt (2 . 0) *s31} . 

Below  I  convert  the  incoming  stresses  to  MY  notation 
Also  note  that  in  the  case  of  use  with  plane-strain  elements 
only  4  components  of  the  stress  are  provided  by  ABAQUS  and  these 
are  the  11,22,33  and  12  components.  In  my  notation,  plane  calculations 
are  performed  in  the  2-3  plane.  So  care  has  to  be  taken 
in  converting  to  my  notation, 
if  (ntens.eq.4)  then 
sc(l)=stress(3) 
sc(2)=stress (1) 
sc(3)=stress{2) 
sc (4) =0 . 0 

sc (5) =dsqrt (2 .DO) *stress (4) 
sc (6) =0 . 0 
else 

sc(l)=stress(l) 
sc(2)=stress(2) 
sc(3)=stress(3) 
sc(4)=dsqrt(2.D0) *stress (4) 
sc { 5 ) =dsqrt ( 2 .DO ) *stress ( 6 ) 
sc { 6 ) =dsqrt ( 2 . DO ) *stress ( 5 ) 
endif 

s(l,l)=sc(l) 

s(2,2)=sc(2) 

s(3,3)=sc(3) 

s(l,2)=sc(4) /dsqrt(2.D0) 
s(2,3)=sc(5) /dsqrt(2.D0) 
s (3, 1) =sc (6) /dsqrt (2 .DO) 
s(2,l)=s(l,2) 
s(3,2)=s(2,3) 
s(l,3)=s(3,l) 

•  Decompose  deltaF  to  get  R  and  U 
Again,  we  must  convert  dfgrdO  to  my  notation: 
if  (ntens.eq.4)  then 
dummyl (1,1) =df grdO (3,3) 
dummyl (2,2) =df grdO (1,1) 
dummyl (3,3) =df grdO (2,2) 
dummyl (1,2) =df grdO (3,1) 
dummyl (2,1) =df grdO (1,3) 
dummyli 2,3) =df grdO (1,2) 
dummyl (3,2) =df grdO (2,1) 
dummyl (1,3) =df grdO (3,2) 
dummyl (3,1) =df grdO (2,3) 
else 

do  10  i=l,3 
do  20  j=l,3 

dummyl ( i , j ) =df grdO ( i , j ) 
continue 
continue 
endif 

call  inverse3x3 (dummyl, dummy2) 
if  (ntens.eq.4)  then 
dummyl ( 1 , 1 ) =df grdl ( 3 , 3 ) 
dummyl (2,2) =df grdl (1,1) 
dummyl (3,3) =df grdl (2,2) 
dummyl {1,2) =df grdl (3,1) 
dummyl (2,1) =df grdl (1,3) 
dummyl ( 2 , 3 ) =df grdl ( 1 , 2 ) 
dummyl (3,2) =df grdl (2,1) 
dummyl (1,3) =df grdl (3,2) 


dummyl (3,1) =df grdl (2,3) 
else 

do  11  i=l,3 
do  21  j=l,3 

dummyl (i, j ) =dfgrdl (i, j ) 

21  continue 

II  continue 
endif 

call  matprod  ( dummyl ,  dummy2 ,  dummy3 ) 
call  decompose (dummy 3, R, strain, eigen, eigval) 
c  write(16,*) 'dfgrdll' ,dfgrdl (1, 1) ,dfgrdl(l,2) ,dfgrdl(l,3) 

c  V7rite(16,*) ' dfgrdl2 ' , dfgrdl (2 , 1) ,dfgrdl(2,2) ,dfgrdl(2,3) 

c  wri te ( 16 , * ) ' df grdl3 ' , dfgrdl (3,1), dfgrdl (3,2), dfgrdl (3,3) 

c  write(16,*) 'strainl' ,strain(l,l) ,strain(l,2) ,strain{l,3) 

c  write (16, *) 'strain2' ,strain(2, 1) , strain (2, 2) , strain (2, 3) 

c  write(16,*) 'strain3 ', strain (3 , 1) ,strain(3,2) ,strain(3,3) 

c _ Here  we  have  obtained  ln(delta_U)  and  R,  the  coirponents  of  which  are 

c  relative  to  the  fixed  Lab  coords, 

c  write (15, *) ' eigen' , eigen 

c  call  flushdS) 

eil (1) =eigen(l, 1) 
eil(2)=eigen(l,2) 
eil (3) =eigen(l, 3) 
ei2(l)=eigen(l,2) 
ei2 (2) =eigen(2, 2) 
ei2  (3) =eigen(3,2) 
ei3  (1) =eigen(l, 3) 
ei3 (2) =eigen(2, 3) 

ei3  (3)  =eigen(3,3)  , 

cwrite(15, *) ' testl=' ,eil (1) *ei2 (1) +eil (2 ) *ei2 (2) +01! (3) *ei2 (3) 
cwrite(15, ♦) ' test2=' , eil (1) *ei3 (1) +eil (2) *ei3 (2) +eil (3) *ei3 (3) 
cwritedS,*)  'test3  =  ' ,ei2  (1)  *ei3  (l)+ei2  (2)  *ei3  (2)+ei2  (3)  *en  (3) 
cwritedS,  *)  '  test4= ' ,  eil  (1)  *eil  (1) +eil  (2)  *eil  (2) +eil  (3)  *eil  (3) 
cwritedS,  *)  '  test5=' ,  ei2  (1)  *ei2  (l)+ei2  (2)  *ei2  (2)  +ei2  (3)  *ei2  (3) 
cwritedS,  *)  '  test6=' ,  ei3  (1)  *ei3  (l)+ei3  (2)  *ei3  (2)  +ei3  (3)  *ei3  (3) 
c  write (15 ,*)' t7 ', eigval (1) , '  ' , eil (1) , eil (2 ) , eil (3) 

c  write (15, *) ' t8' , eigval (2) , '  ' ,ei2 (1) , ei2 (2) , ei2 (3) 

c  write (15, *) ' t9 ' , eigval (3) , '  ' , ei3 (1) , ei3 (2) , ei3 (3) 

if  ( (kinc.eq.l) .and. (kstep.eq.l) )  then 
f=0.15 
wil=l . 0001 
wi2=l . 0001 
atheta=0 . 0 
aphi=0 . 0 
apsi=:0 . 0 
tottheta=0.0 
totphi=0.0 
totpsi=0 . 0 
do  777  i=l,3 
do  778  j=l,3 
if  (i.ne.j)  then 
rotj (i, j ) =0 .0 
else 

rotj (i, j)=1.0 
endif 

778  continue 

III  continue 
statevd)  =:f 
sigyl=props (5) 
statev(2) =wil 
statev(3) =wi2 
statev ( 4 ) =atheta 
statev(5) =aphi 
statBv(6)  =:apsi 
statev(7)=0.0 
statev ( 8 ) =props ( 5 ) 

statev (9) =0 . 0 
statev  dO)  =hrd 
‘ statev (11) =wil/wi2 
statev ( 12 ) =rot j (1,1) 
statev (13 ) =rot j (1,2) 
statev ( 14 ) =rotj (1,3) 
statev (15 ) =rotj (2,1) 
statev (16 ) =rotj (2,2) 
statev (17 )=rotj (2,3) 
statev ( 18 ) =rotj (3,1) 
statev (19 ) =rotj (3,2) 
statev ( 20 )=rotj (3,3) 
statev (21 ) =tottheta 
statev (22 ) =totphi 
statev (23 ) =totpsi 
statev (24) =1 . 0/wi2 
statev (25) = (1 . 0/wi2) *wil 
endif 
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do  5  i=l,3 
do  6  j=l,3 

if  (dabs(R(i, j) ) .lt.l.Od-14)  then 
R(i, j)=0.0 
endif 
continue 
continue 
do  100  i=lt3 
do  110  j=l,3 
RT(i, j)=R(j,i) 
continue 
continue 


-  This  is  a  test  to  see  if  the  components  of  [drot] [s]_t tdrot_transpose] 
that  ABAQUS  gives  me  are  with  respect  to  the  fixed  frame.  So,  what  I 
do  below  is  assume  that  this  is  true  (they  are  with  respect  to  the 
GLOBAL  axes)  and  then  use  transformation  laws  to  obtain  the  components 
w.r.t  to  a  frame  which  is  obatined  from  the  GLOBAL  frame  by  a  rotation 
of  drot .  These  con^onents  are  the  same  as  the  components  of  [s]_t  w.r.t 
the  GLOBAL  frame.  I  then  check  if  these  values  agree  with  the  stress 

of  the  previous  increment:  this  is  output  in  the  .dat  file  and  since 

all  tensors  are  stored  w.r.t  the  GLOBAL  axes,  the  stress  output  should 

agree  with  'dummy2'  below  and  they  do! 

write (15, *) 'stress^l' 

duinmy2  (1, 1)  =stress  (1) 

dummy2 (2, 2) =stress (2) 

d\iinmy2  (3,3)  =stress  ( 3 ) 

durnmy2  (1,2)  =s tress  (4) 

dummy2 (2,1) =s tress (4) 

dummy2 (2,3) =stress (6) 

dummy2 (3,2) =stress (6) 

dummy2 (3,1) =stress ( 5 ) 

duinmy2  ( 1 ,  3 )  =stress  ( 5 ) 

call  inatprod(RT,dummy2,dummyl) 

call  matprod(dummyl,R,dummy2) 

write  (15,  *)  duinmy2 

-  END  TESTl 

-  BEGIN  TEST2 

-  From  dfgrdO  and  dfgrdl,  I  have  obtained  ln(delta_JJ)  and  R.  R  is  exactly 
the  same  as  drot.  The  dstran  that  ABAQUS  provides  me  is  ln(delta_V)  and 
its  components  are  w.r.t  the  GLOBAL  axes.  I  have  the  components  of 
ln(delta_U)  also  w.r.t  the  GLOBAL  axes.  It  can  be  easily  checked  to 
see  that  they  have  the  same  eigen  values.  Also  [R] [In (delta_U) ] [RT] 
must  give  us  dstran. 

writedS,  *)  'R' 
write (15, *)R 
write(15,*) 'drot' 
write (15, *)drot 
do  31  i=l,3 
do  32  j=l,3 

dummyl ( i , j ) =strain ( i , j ) 
continue 
continue 
n=3 
np=3 

call  jacobi (dummyl, n, np, eigval,duinmy2,nrot) 
write (15,  *) 'Eigen  values  of  ln(delta__U)  ' 
write(15, *)eigval 
do  33  i=l,3 
do  34  j=l,3 

dummyl ( i , j ) =s train { i , j ) 
continue 
continue 

write (15, *) 'ln(delta_U) ' 

write (15,  *) dummyl  (1, 1)  , dumiry  1  ( 2 , 2 )  ,dummyl  (3, 3)  ,dummyl  (1,2)  *2.0, 

&dummyl  (3 , 1)  *2 . 0, dummyl  (2, 3)  *2 . 0 
call  matprod{RT, duranyl,dummy2) 
call  matprod(diHnmy2,R,diiinmyl) 
write(15, *) 'ln(delta_V)  from  ln(delta_U) ' 

write(15, *) dummyl (1,1) , dumn^l ( 2 , 2 ) ,dummyl(3,3) , dummyl ( 1 , 2 ) *2.0, 

&dummyl  (3, 1)  *2 .0,duinmyl  (2, 3)  *2 . 0 
write (15, *) 'ln(delta_V) ' 
write (15, *) dstran 
--  END  TEST2 

-  Stress  conponents  that  abaqus  provides  me  are  w.r.t  GLOBAL 
co-ordinate  frame.  Since  this  is  a  finite -deformation  problem, 

the  stress  that  I  have  here  is  one  that  has  been  rotated  by  'drot', 
i.e,  it  is  [drot]  [stress]_t  [drot_transpose] ,  but  the  components 
are  w.r.t  the  fixed  GLOBAL  frame.  The  integration  algorithm  that 
I  use  requires  that  I  used  [stress]_t  euid  NOT  the  above  quantity. 

This  is  done  below: 

.call  matprod(RT,s, dummyl) 
call  matprod (dummyl , R, s) 


c 

C' 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c- 

c 

c 

c 

c 

c 

c 

c 

99 


— Note  that  R  and  drot  are  exactly  the  same.  The  components  of  s  are 
still  relative  to  the  global  frame. 

An  important  interpretation  of  why  we  use  the  stress  in  the  GLOBAL 
coordinate  frame  is  that  when  we  add  to  this  the  change  in  stress 
we  get  the  new  stress  components  relative  to  a  coordinate  frame 
which  is  obtained  from  the  global  frame  by  a  rotation  of  drot .  This 
stress  is  called  sig_hat_n+l  in  Dr.  Aravas's  notes.  We  see  that  to 
recover  stress_n+l,  we  must  pre  and  post  multiply  by  R  and  RT, 
respectively.  This  is  nothing  but  expressing  the  new  stress  components 
relative  the  GLOBAL  frame. 

Note  that  we  need  not  use  the  GLOBAL  frame  as  one  where  all  components 
are  stored.  In  fact,  we  find  it  convenient  to  perform  our  integration 
relative  to  a  coordinate  frame  which  coincides  with  the  orientation 
of  the  particles.  What  this  implies  is  that,  we  transform  the  stress 
we  have  obtained  above  to  this  (particle  orientation  frame)  coordinate 
frame  and  perform  our  integration  there.  The  stress  components  that  we 
then  obtain  are  then  w.r.t  a  coordinate  frame  which  has  is  obtained 
from  the  particle  orientation  frame  by  a  rotation  of  drot.  We  then 
again  transform  the  stress  (and  other  variables)  to  be  expressed 
relative  to  the  GLOBAl  frame. 

sc(l)=s(l,l) 

sc(2)=s(2,2) 

sc(3)=s(3,3) 

sc (4) =dsqrt (2 .DO) *s (1, 2) 
sc ( 5 ) =dsqrt (2.D0)*s(2,3) 
sc (6) =dsqrt (2.D0)*s(3,l) 
if  ( (kinc.gt.l) .or. (kstep.gt.l) )  then 
f =statev (1) 
wil=statev(2) 
wi2=statev (3) 
c=1.000 
a=c/wil 
b=c/wi2 

atheta=statev (4) 
aphi=statev(5) 
apsi=statev(6) 
eps_plastic=statev(7 ) 
sigyl=statev (8) 
rotj (1, 1) =statev(12) 
rotj (1,2) =statev(13) 
rotj (1,3) =statev(14) 
rotj (2, 1) =statev(15) 
rotj (2,2) =statev(16) 
rotj (2,3) =statev(17) 
rotj (3,l)=statev{18) 
rotj (3,2)=statev(19) 
rotj {3,3)=statev(20) 
tottheta=statev(21) *3 , 14159/180 . 0 
totphi=statev (22) *3 . 14159/180 . 0 
totpsi=statev(23) *3.14159/180.0 
if  (f. It. 0.001)  then 
evolf=, false. 
evolwi= . false . 
endi  f 

if  (wi2.lt. 0.02)  then 
evolwi= . false . 

endif 

endif 

--  The  increment  of  strain  below  (corresponds  to  ln(delta_U))  is  w.r.t 
GLOBAL  coordinate  frame, 
strainc ( 1 ) =s train (1,1) 
strainc (2) =strain(2, 2) 
strainc ( 3 ) =strain (3,3) 
strainc (4) =strain (1, 2) *dsqrt (2 .DO) 
strainc (5) =strain (2 , 3) *dsqrt (2 .DO) 
strainc ( 6 ) =strain (3,1) *dsqrt ( 2 . DO ) 

Expressing  stress  euid  strain  in  local  coordinates: 

The  integration  problem  is  carried  out  in  a  co-ordinate  frame 
which  coincides  with  the  orientation  of  the  particles.  This  is 
done  because  it  is  more  eonomical  to  rotate  the  stress  and  the 
strain  once  rather  than  having  to  rotate  4th  order  tensors  like 
Amat,  Bmat,  MHS  etc.  whose  components  are  readily  obtained  w. r . t ^ 
the  co-ordinate  frame  which  coincides  with  the  particle  orientation. 

ct=dcos (atheta) 
st=dsin(atheta) 
cp=dcos (aphi) 

cs=dcos (apsi)  ' 

sp=dsin(aphi) 
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ss=dsin (apsi) 
rota  (1, 1)  =cs*cp-'ss*ct*sp 
rota (1 , 2 ) =-cs*sp-ss*ct*cp 
rota(l,3)=ss*st 
rota (2 , 1) =ss*cp+cs*ct*sp 
rota (2,2) =-ss*sp+cs*ct*cp 
rota(2,3)=-cs*st 
rota (3 , 1) fst*sp 
rota(3,2)=st*cp 
rota(3,3)=ct 

call  matprod (rota, rot j , rot) 
do  7  i=l,3 
do  8  j=l,3 
rott(i, j)=rot(j,i) 
continue 
continue 

-  Below  the  components  of  the  stress  are  expressed  w.r.t.  the 
a  coordinate  system  which  instantaneously  coincides  with  the 

c  orientation  of  the  particles. 

call  inatprod(rott,  s,dvimmyl) 

call  ma tpr od( dummy 1, rot, s) 

sc(l)=s(l,l) 

sc(2)=s(2,2) 

sc(3)=s(3,3) 

sc (4) =dsqrt (2 .DO) *s (1, 2) 
sc (5) =dsqrt (2 .DO) *s (2, 3) 
sc { 6 ) =dsqrt { 2 . DO ) *s ( 3 , 1 ) 

c _  Rotating  below  to  express  strain  w.r.t  the  orientation  of  particles. 

call  matprod(rott, strain, dummyl) 
call  ma  tpr  od(diimmyl,  rot,  strain) 
strainc (1) =strain(l, 1) 
strainc (2) =s train (2,2) 
strainc (3) =strain(3,3) 
strainc (4) =s train (1, 2) *dsqrt (2 .DO) 
strainc (5) =strain(2 ,3) *dsqrt (2 .DO) 
strainc ( 6 ) =strain (3,1) *dsqrt ( 2 . DO ) 

Q -  Integrating  to  obtain  the  stress  and  the  new  state  variables. 

xacc=l . Od-9 

if  {  (strained)  .eq. 0.0)  .and.  (strainc(2)  .eq. 0.0) 

&  .and. (strainc(3) .eq.0.0) .and. (strainc(4) .eq.0.0) 

Sc  .and.  (strainc(5)  .eq.0.0)  .and.  (strainc(6)  .eq.0.0) )  then 
dlam=0.0 
go  to  4 
endif 

call  guessmaker (x2) 
if  (yield. eq. . false. )  then 
dlam=0 . 0 
go  to  4 
endif 

if  (neg. eq. . true . )  then 
pnewdt=0 .75 

write (15, *) 'This  should  never  happen!' 
call  flushdS) 
return 
endif 

if  {x2.lt. 0.0)  then 

c  wr ite(l 5, *) 'Guessmaker  screwed  up',x2 

c  call  flush(15) 

endif 

Q -  Trying  to  bracket  a  solution  for  delta_lambda 

x2=0.02 
xl=1.0d-10 
num=20 
nb=l  ■ 

call  zbrak(func,xl,x2,num,xbl,xb2,nb) 
if  (neg. eq. .true. )  then 
pnewdt=0 .75 
return 
endif 

if  (nb.eq.O)  then 

c  wri te(l 5, *) 'bracketing  problem:  Level  1' 

call  flush (15) 
xl=1.0d-10 
num=50 
nb=l 

call  zbrak ( f unc , xl , x2 , num, xbl , xb2 , nb) 
if  (neg. eq. .true. )  then 
pnewdt=0.75 
return 
endi  f 
endi  f 

if  (nb.eq.O)  then 

write (15, *) 'bracketing  problem:  Level  2' 
call  flushdS) 


c 


zb= . true . 
xl=1.0d-10 
n\im=200 
nb=l 

call  zbrak ( f unc , xl , x2 , num, xbl , xb2 , nb) 
if  (neg.eq. . true. )  then 
pnewdt=0 .75 
return 
endif 
endif 

if  (nb.eq.O)  then  ,  , 

write (15, *) 'bracketing  problem:  Level  critical' 
call  flush (15) 
xl=1.0d-10 
num=500 
nb=l 

call  zbrak (func,xl,x2, num, xbl, xb2,nb) 
if  (neg.eq. . true. )  then 
pnewdt=0 . 75 
return 
endif 
endif 
zb=. false. 

if  (nb.eq.O)  then 
pnewdt=0 .75 

write (15, *) 'Bracketing  Failure' 
call  flush(15) 


return 

c _  Using  a  combined  secant/bisection  procedure  to  find  the  solution  for 

c  delta_lambda. 

c 

dlam=rtsaf e ( f uncd , xbl , xb2 , xacc) 
c  write (15, *) 'dlam' ,dlam, 'kinc' ,kinc 

if  (neg.eq. .true. )  then 

pnewdt=0,75  n 

write (15, *) 'returned  after  unsuccessful  rtsafe' ,kinc,noel,npt 

call  flush(15) 


return 

endif 

c _  Initializing  theta;  The  voids  start  evolving  at  an  orientation  given 

c  by  the  eigen  values  of  the  right-stretch  tensor.  The  orientation  of  the 

c  voids  makes  sense  only  when' the  voids  are  not  spherical  anymore.  This 

c  means  that  the  following  steps  are  carried  out  only  when  the  material 

c  has  started  deforming  plastically, 

if  (evolti.eq. .true. )  then 
if  (statevO)  .eq.0. 0)  then 
write (15, *) 'ei2 ' , ei2 
write(15, *) 'ei3' ,ei3 
cwrite(15, ♦) 'thetaj ' ,dacos(rotj (3,3) ) 
call  flush(15) 

if  (dabs(ei2(3) ) .gt.l.Od-10)  then 
if  (eigval (3) .gt . eigval (2) )  then 

if  ( (ei2(3) .lt.0.0) .and. (ei2(2) .gt.0.0) )  then 
atheta=-dacos (ei2 (2) ) 

write (15, *)  ' theta_initializing_a' ,noel,npt 

else  if  ((ei2(3) .gt.0.0) .and. (ei2(2) .gt.0.0))  then 
atheta=dacos (ei2  (2) ) 

write (15, *)  ' theta_initializing_b' ,noel,npt 

else  if  ( (ei2(3) .gt.0.0) .and. (ei2(2) .lt.0.0) )  then 
write ( 15 , * )  ' theta_ini tializing : warningl ' 

atheta=dacos (ei2  (2) ) 

else  if  ( (ei2(3) .lt.0.0) .and. (ei2(2) .lt.0.0) )  then 
artheta=-dacos  (ei2  (2)  ) 
write ( 15 , * )  ' theta_initializing : warning2 ' 

else 

atheta=0.0 

endif 

else  if  (eigval (3) . It. eigval (2) )  then 
if  ((ei2(3) .lt.0.0) .and. (ei2(2) .gt.0.0))  then 
atheta= ( 3 . 14159/2 .0) -dacos (ei2 (2 ) ) 
write (15, *)  ' theta_initializing_c' ,noel,npt 
else  if  ( (ei2(3) .gt.0.0) .and. (ei2 (2) .gt.0.0) )  then 
atheta=- ( (3 . 14159/2 . 0) -dacos (ei2 (2) ) ) 
write (15, *)  ' theta_initiali zinged' ,noel,npt 

else  if  ( (ei2 (3) .gt.0.0) .and. (ei2 (2) .lt.0.0) )  then 
write ( 15 , * )  ' theta_initializing : warning3 ' 

atheta=-dacos ( ei2 ( 2 ) ) 

else  if  ( (ei2(3) .lt.0.0) .and. (ei2 (2) .lt.0.0) )  then 
atheta=dacos (ei2 (2) ) 

write (15 , * )  ' theta^initializing ; warning4 ' 

else 

atheta=0 .0 
endif 


u 


else 

write (15, *) 'voids  are  spherical : problems ' 
endif 

atheta=0 . 0 

statev  ( 9 )  =statev  { 9*)  +1 . 0 
write (15, *) 'sending  back' 
call  flush (15) 
go  to  99 
else 

atheta=0 . 0 

statev ( 9 ) =statev ( 9 ) +1 . 0 
endif 
endif 

^Using  the  solution  for  'dlaiti'  (delta_lainbda)  to  update  all  variables 


call  update (dlam,ddsdde_mine) 

if  ( (f .lt.0.0) .or . (wil.lt. 0.0) .or. (wi2.lt. 0.0) )  then 

pnewdt=0.75  .  u  -u 

write (15, *) 'This  should  never  happen:  problem  should  have  been 
Sc  taken  care  of  before  updating' , dlam, f ,wil,wi2 
call  flush (15) 
return 
endif 

_  stress  and  ddsdde  are  relative  to  the  GLOBAL  axes . 

s(l,l)=sc(l) 

s(2,2)=sc(2) 

s(3,3)=sc(3) 

s(l,2)=sc(4) /dsqrt(2.D0) 

s(2,3)=sc(5) /dsqrt(2.D0) 

s(3,l)=sc(6) /dsqrt(2,D0) 

s(2,l)=s(l,2) 

s(3,2)=s(2,3) 

s(l,3)=s(3,l) 

sc(l)=s(l,l) 

sc(2)=s(2,2) 

sc(3)=s(3,3) 

sc (4) =s (1,2) *dsqrt(2.D0) 
sc(5)=s(2,3)*dsqrt(2.D0) 
sc(6) =s(3, 1) *dsqrt (2 .DO) 
if  (ntens.eq.4)  then 
stress (1) =sc  (2) 
stress (2 ) =sc (3 ) 
stress (3) =sc (1) 
stress (4) =sc (5) /dsqrt (2 .DO) 
else 

stress (1) =sc (1) 
stress (2) =sc  (2) 
stress (3) =sc (3) 
stress (4) =sc (4) /dsqrt (2 .DO) 
stress (5) =sc (6) /dsqrt (2 .DO) 
stress (6) =sc (5) /dsqrt (2 .DO) 
endif 

statev(l)=f 

if  ( (statev(3) .gt.1.0) .and. (tottheta. It . 0 . 0) )  then 
wi2=l . 0/wi2 
wil=wi2*wil 

tottheta=0 . 5*3 . 14159+tottheta 
endif 

statev  (2)  =wil 

statev  (3 )  =wi2 

s  tatev ( 4 ) =athe ta 

statev (5) =aphi 

statev(6) =apsi 

s tatev(  7 ) =eps_plastic 

statev ( 8 ) ^sigyl 

statev ( 10 ) =hrd 

statev (11) =wil/wi2 

statev ( 12 )=rotj (1,1) 

statev ( 13 ) =rotj (1,2) 

statev (14 ) =rotj (1,3) 

statev ( 15 ) =rotj (2,1) 

s tatev ( 16 ) =rotj (2,2) 

statev(17) =rotj (2,3) 

statev (18 ) =rotj (3,1) 

statev (19 ) =rotj (3,2) 

statev (20) =rotj (3,3) 

tot theta=dasin ( rot (3,2)) 

statev (21 ) =tottheta*180 .0/3 .14159 

statev (22 ) =totphi* 180 .0/3.14159 

statev (23 )=totpsi*180. 0/3. 14159 

statev (24 ) =1 . 0/wi2 

statev (25)  =  (1 .0/wi2)  *wil 

statev(26) =1 . 000/wil 

-  converting  ddsdde  for  the  case  of  2D  analyses: 


if  (ntens.eq.4)  then 
ddsdde  ( 1  #  1 )  =ddsdde_itiine  (2,2) 
ddsdde (1,2) =ddsdde_mine (2,3) 
ddsdde  (1,3)  =ddsdde_inine  (2,1) 
ddsdde  (1,4)  =ddsdde__niine  (2,6) 
ddsdde (2,1) =ddsdde_mine (3,2) 
ddsdde  (2,2)  =ddsdde_iiiine  (3,3) 
ddsdde (2,3) =ddsdde_mine (3,1) 
ddsdde  (2,4)  =ddsdde_inine  (3,6) 
ddsdde  (3,1)  =ddsdde_niine  (1,2) 
ddsdde  (3,2)  =ddsdde_inine  (1,3) 
ddsdde  (3,3)  =ddsdde_inine  (1,1) 
ddsdde  (3,4)  =ddsdde__niine  (1,6) 
ddsdde (4,1) =ddsdde_mine (6,2) 
ddsdde (4,2) =ddsdde_mine (6,3) 
ddsdde (4,3) =ddsdde_mine (6,1) 
ddsdde  ( 4  /  4 )  =ddsdde_inine  (6,6) 
else 

do  i=l,6 
do  j =1,6 

ddsdde ( i , j ) =ddsdde_mine ( i , j ) 
enddo 
enddo 
endif 
111  return 
end 


subroutine  funcd(dlaiti,  fun,  dfun) 
real *8  fun, func, dlam, dfun 
logical  debug, neg, check, yield 
common  /mkdata2 /yield, check, neg 
common  /mdebug/ debug 
external  func 


fun=func (dlam) 

if  (neg. eq. .true. )  then 

write (15, *) 'Returning  from  FUNCD:  Problem  due  to  func  evaluation' 
call  flushdS) 
return 
endif 

dfun= ( func (dlam+0 . 0001*dlam) -fun) / (0 . 0001*dlam) 
if  (dfun. eq, 0 . 0)  then 

write (15, *) 'You  should  never  see  this  statement:  Error  trapping 
Sc  incorrect  in  func' 
call  flushdS) 
neg= . true . 
endif 

if  (neg. eq. .true. )  then 

write (15, ♦) 'Returning  from  FUNCD.  Problem:  In  func  during 
Sc  slope  calculation' 
call  flushdS) 
return 
endif 

if (debug .eq. . true. )  then 

write (15, * ) 'dfun' , dfun, fxin, func (dlam+0 . 0001*dlam) 

endif 

return 

end 


subroutine  update (dlam, ddsdde) 
real*8  s(3,3) ,MHS(6,6) ,Amat(6,6) ,Ce(6,6) 
real*8  sn(3,3),n(3,3) , omega (3, 3) ,sige(3,3) ,L 
real * 8  strainc ( 6 ) , sigec ( 6 ) , ddsdde (6,6) , snl (3,3) 
real*8  sigma_n(6) , sig_n(6) ,ddsdum(6, 6) 

■  real *  8  f , wil , wi2 , athe ta , aphi , aps i , athetan , aphin , aps in 
real *8  fn,  wiln,  wi2n,  dlam,  durnm(3 , 3,3,3) 
real*8  sc (6) , nc (6) , omegac (6) , sigyln 
real*8  mul,kl,mu2,k2, sigyl,H,dumml (3, 3, 3, 3) 
real*8  dummy(6) , a,b, c, ad, bd, cd, dummyl (3 , 3 ) ,dummy2(3,3) 
real *8  dummyS (6) , duml ,nul,nu2 , sigyO 
real * 8  deps (3,3), eps_plas tic , dep_j)las , dumd (6,6) 
real*8  det,kll,dphidf ,phil,phi2, ftest 

real*8  el212,e2323,el313,pil212,pi2323,pil313,bmat (6, 6) 

real* 8  mule, kle,mu2e, ek2e, Ce_temp (6, 6) 

real * 8  ate , bte , dphidwil , dphidwi2 , wil test , wi2 test 

real *8  dum5,duml0,y,dinc (6) ,R(3,3) ,RT(3,3) 

real*8  cp, sp, ct, st, cs, ss, rot (3, 3 ) ,rott(3,3) 

real*8  cpn, spn, ctn, stn, csn, ssn, rotn (3 , 3 ) , rotnt (3 , 3 ) 

real*8  vec(3) ,hrd,rotj (3,3) ,rotjn(3,3) ,rota(3,3) 

real*8  test (6 , 6) , loadp (6) , loading 


real *8  zero ( 6) , tots train (6) 

integer  i, j ,ninc, indx(6) , inpt 

logical  path, yield, check, neg, spath,yc, debug 

logical  evolf , evolwi, evolti 

common  /controls/ evolf, evolwi, evolti 

common  /mJonodulus  /mul ,  kl ,  mu2 ,  k2 ,  s  igy 0 ,  kl  1 

common  /mkelas/mule,  kle,mu2e,  ek2e 

common  /mkdatal/f  ,wil, wi2,  sc,  strainc, eps_plastic,  sigyl,hrd  ^ 

common  /mkorient/cp, sp, ct , st , cs , ss, rot, rott , atheta, aphi , apsi 
common  /mkdata2 /yield, check, neg 
common  /mkrot/R,RT, rotj 
common  /paths/ spa th 
common  /mdebug/ debug 
c 

path= . true . 
neg=. false, 
c 

i f ( debug . eq . . true . )  then 
write (15, *) 'updatel' 
call  flush (15) 
endif 

c  write ( 16, *) 'strainc' , strainc 

nul=0.5* (3 .0*kl-2.0*mul) / (3 .0*kl+mul) 
nu2=0.5* (3.0*k2-2.0*mu2) /(3.0*k2+mu2) 
c=1.000 
a=c/wil 
b=c/wi2 
ad=a 
bd=b 
cd=c 

if  (yield. eq. .false. )  then 

c -  The  elastic  predictor  is  the  correct  stress 

c  The  state  variables  are  unchanged  by  elastic  deformations 

c  and  so  there  is  no  need  to  update  them. 

path=. false. 

call  Mef fective (a,b, c, ad,bd, cd, f, mule, kle,mu2e, ek2e, Ce, path) 
path= . true . 

call  teninatprod(Ce, strainc, sigec) 
sige (1,1) =sc ( 1 ) +sigec ( 1 ) 
sige(2,2) =sc(2) +sigec(2) 
sige (3 , 3 ) =sc  O) +sigec (3) 
sige(l,2)=(sc(4)+sigec(4) ) /dsqrt (2.D0) 
sige (2,3) = (sc (5) +sigec (5) ) /dsqrt (2 .DO) 
sige(3,l)=(sc(6)+sigec(6) ) /dsqrt (2 .DO) 
sige (2 , 1) =sige (1, 2) 
sige (3 , 2) =sige (2 , 3) 
sige (1, 3 ) =sige (3 , 1) 

c -  update  stress  (state  variables  are  unchanged!) 

Q -  Express  stress  in  GLOBAL  coordinate  frame. 

call  matprod (rot, sige, dummyl) 
call  matprod (dummyl, rott, sige) 
sc(l) =sige(l, 1) 
sc(2)=sige(2,2) 
sc{3)=sige(3,3) 
sc ( 4 ) =sige (1,2) *dsqrt ( 2 . DO ) 
sc (5) =sige (2 , 3) *dsqrt (2 .DO) 
sc(6) =sige{3, 1) *dsqrt (2.D0) 


c -  calculate  ddsdde 

do  431  i=l,6 
do  441  j=l,6 
dumd ( i , j ) =Ce ( i ,  j ) 
441  continue 

431  continue 


c -  Express  ddsdde  in  GLOBAL  coordinate  frame 

call  mat2tensor  (duind,duinm) 
call  rot4order  (dumm,  rot,duinml) 
call  ten2matrix{dumml,Ce_temp) 

c -  Expressing  ddsdde  in  the  notation  of  ABAQUS 

do  432  i=l,3 

dumd(i,  4)  =Ce_temp ( i ,  4 )  /dsqrt  (2  .DO) 
dumd{i,  5)=Ce_temp(i,5)/dsqrt(2.D0) 
dumd  ( i ,  6 )  =:Ce_teitp  ( i ,  6 )  / dsqrt  ( 2 .  DO ) 

432  continue 

do  433  i=4,6 

dumd ( i , 1 ) =Ce_temp ( i , 1 ) / dsqrt ( 2 . DO ) 
dumd(i,  2)=Ce_teit^(i,2)  /dsqrt  (2.D0) 
dumd ( i , 3 ) =Ce_temp ( i , 3 ) /dsqrt ( 2 , DO ) 

433  continue 

do  434  i=4,6 
do  435  j=4,6 

dumd  ( i ,  j )  =Ce_teii^  ( i ,  j )  /  2 . 0 
continue 
continue 
do  436  i=l,4 


435 

434 


do  437  j=l,4 
ddsdde  {i,  j )  =diimd (i ,  j  ) 
continue 
continue 
do  438  i=l,4 
ddsdde(i,  5)  =duind(i,  6) 
ddsdde U,  6)  =duind(i,  5) 
continue 
do  439  i=i,4 
ddsdde  ( 5 ,  i )  =duind  ( 6 ,  i ) 
ddsdde  ( 6 ,  i )  =duind  ( 5 ,  i ) 
continue 

ddsdde{5,5)=dumd(6,6) 
ddsdde  (5,6)  =duind  (6,5) 
ddsdde (6,5) =dumd (5,6) 
ddsdde (6,6) =dumd (5,5) 
call  ludcmp(Ce_teinp,  6,6,indx,det) 
do  544  j=l,6 
det=det*Ce_temp ( j , j ) 
continue 

if  (det.lt. 0.0)  then 
write (15 , * ) ' ddsdde_elastic ' , a, b, c 
write (15, *) ddsdde 
write(15,*) 'det',det 
endif 
go  to  999 

endif  .  . 

If  the  material  has  become  plastic,  the  microstructural  variables 
need  to  be  updated.  Also  the  new  stress  has  to  be  evaluated  carefully 
since  the  strains  are  no  longer  completely  elastic.  We  must  obtain 
the  plastic  part  of  the  strain  (which  is  done  by  obtaining  'dlam' 
i.e.,  delta^lambda)  and  then  that  is  used  to  update  all  relevant 
variables . 

-  calculate  Meffective 

call  Meffective(a,b,c,ad,bd,cd, f ,mul,kll,mu2,k2,MHS,path) 
do  100  i=l,6 
do  110  j=l,6 

MHS(i, j)=3.0*mul*MHS(i, j) 
continue 
continue 


estimate  stress 

--  STEP  1:  Estimate  the  elastic  predictor 
path=. false. 

call  Meffective (a, b,c, ad, bd, cd, f ,mule, kle,mu2e, ek2e,Ce,path) 
path= . true . 

call  teninatprod(Ce,  strainc,  sigec) 
do  201,  i=l,6 
loadp(i)=sigec(i) 
continue 

sige (1,1) =sc ( 1 ) +sigec ( 1 ) 

sige (2 , 2 ) =sc (2) +sigec (2) 

sige (3, 3 ) =sc (3) +sigec (3) 

sige (1 , 2 ) = (sc (4 ) +sigec (4) ) /dsqrt (2 .DO) 

sige (2, 3 ) = (sc (5) +sigec (5) ) /dsqrt (2. DO) 

s ige ( 3 , 1 ) = ( sc ( 6 ) +s igec ( 6 ) ) /dsqrt ( 2 . DO ) 

sige (2, 1) =sige(l,2) 

sige (3,2) =sige (2 , 3 ) 

sige ( 1 / 3 ) =sige (3,1) 

sigec(l) =sige(l,  1) 

sigec (2) =sige (2 , 2) 

sigec(3) =sige(3, 3) 

sigec ( 4 ) =sige (1,2) *dsqrt ( 2 . DO ) 

sigec(5) =sige (2,3) *dsqrt (2 .DO) 

sigec (6) =sige (3, 1) *dsqrt (2 .DO) 

—  STEP  la:  Determining  the  stress  on  the  yield  surface 
sigma_n ( 1 ) =sc ( 1 ) 
sigma_n(2) =sc (2) 
sigma_n ( 3 ) =sc ( 3 ) 
sigma_n  (4)  =sc  (4) 
sigma_n(5) =sc(5) 
sigma_n ( 6 ) =sc ( 6 ) 

call  yieldtes t  (MHS ,  sigma_n ,  sigyl ,  f ,  yc ,  y ) 
if  (yc.eq. . false. )  then 
call  yieldtest (MHS , sigec , sigyl , f i yc , y ) 
if  (yc.eq. .true. )  then 

this  is  the  step  where  the  behavior  becomes  plastic  (the  previous 
step  was  elastic) 
check= . true . 
do  127  i=l,6 
zero (i) =0 . 0 

totstrain ( i ) =strainc ( i ) 
continue 
ninc=100 . 0 
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c 


203 


c 

c 


42 

41 


call  stre (sigina_n, sigec, zero, totstrain, sig_n, Ce, nine) 
else 

write{15,  *) 'NOT  PLASTIC,  shouldnt  have  come  here 
call  flush (15) 
do  130  i=l,6 
sig_n{i)=sc(i) 
continue 
endif 
else 

do  131  i=l,6 
sig_n(i)=sc(i) 
continue 
endif 

—  calculate  N(3,3) 

call  tenmatprod(MHS,sig_n,nc) 
do  200  i=l,6 
nc(i) =nc (i) / (1. 0-f ) 
nc(i)=2.0*nc(i) /sigyl 
dummy ( i ) =nc ( i ) 
dummy 5 ( i ) =nc ( i ) 
continue 

write ( 16 , * ) 'mhs ' , mhs (1,1) , mhs (1,2) , mhs (1,3) , mhs (1,4) 

Sc  ,mhs  (1, 5)  ,inhs  (1,  6) 

write (16,*) 'mhs2 ' ,mhs (2 , 1) , mhs (2, 2) ,mhs(2,3) ,mhs(2,4) 

Sc  ,mhs  (2, 5)  ,mhs  (2, 6) 

write (16, *) 'sign' ,sig_n 
write (16, *) 'nc' ,nc 
write (15, *) 'exit  stress' 
write (15 , *) sig_n 
n(l,l)=nc(l) 
n(2,2)=nc(2) 
n(3,3)=nc(3) 

n(l,2)=nc(4) /dsqrt(2.D0) 
n(2,l)=n(l,2) 
n(2,3)=nc(5) /dsqrt(2.D0) 
n(3,2)=n(2,3) 
n(3,l)=nc(6) /dsqrt(2.D0) 
n(l,3)=n(3,l) 

—  Test  for  loading/unloading 

call  rowcolumnprod(loadp, dummy5, loading) 
do  203  i=l,6 
dummy5 ( i ) =nc ( i ) 
continue 

if  ( loading. le . 0 . 0)  then 
write (15, *) ' loading' , loading 
endif 

' —  STEP  2:  Remaining  part  of  the  stress 
- —  STEP  2a:  Calculating  'omega' 

spath=.  false.  . 

call  stensor (a,b,c,nul,kl,mul,dumd,pil212,pil313,pi2323) 

if  ( debug. eq. , true. )  then 

write (15, *) 'pil212 ' ,pil212 ,pi2323 ,pil313 

call  flush(15) 

endif 

el212=pil212-f*pil212 
e2323=pi2323-f*pi2323 
el313=pil313-f*pil313 
spath= . true . 

call  Btensor (a, b, c, ad, bd, cd, f ,mul, kl ,mu2 , k2 , Bmat ) 
do  41  i=l,3 
do  42  j=l,6 
Bmat(i, j)=0.0 
continue 
continue 

Bmat (4,l)=2.0*el212  *Bmat (4,1) 
Bmat(4,2)=2.0*el212*Bmat(4,2) 

Bmat (4, 3) =2 .0*el212*Bmat (4,3) 

Bmat (4, 4) =2 . 0*el212*Bmat (4, 4) 

Bmat (4, 5) =2 .0*el212*Bmat (4, 5) 

Bmat (4, 6)=2.0*el212*Bmat(4,6) 

Bmat ( 5 , 1) =2 . 0*e2323*Bmat ( 5, 1) 

Bmat (5, 2) =2 .0*e2323*Bmat (5, 2) 

Bmat ( 5 , 3 ) =2 . 0*e2323  *Bmat (5,3) 

Bmat ( 5 , 4 ) =2 . 0*e2323*Bmat ( 5 , 4 ) 

Bmat (5, 5) =2 . 0*e2323*Bmat (5, 5) 
Bmat(5,6)=2.0*e2323*Bmat(5,6) 

Bmat (6, 1) =2 .0*el313*Bmat (6, 1) 
Bmat(6,2)=2.0*el313*Bmat(6,2) 
Bmat(6,3)=2.0*el313*Bmat(6,3) 

Bmat ( 6 , 4 ) =2 . 0 *el3 13 *Bmat ( 6 , 4 ) 

Bmat (6,5)=2.0*el313  *Bmat (6,5) 
Bmat(6,6)=2.0*el313*Bmat(6, 6) 

call  Atensor (a,b, c,ad,bd, cd, f ,mul, kl,mu2, k2, Amat) 
do  3  i=l,3 
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do  4  j=l,3 
omega (i, j ) =0 . 0 
4  continue 

3  continue 

do  1  i=l,3 

omega(l,2)=-Binat(4,i) *nc(i)+omega(l,2) 

omega  ( 2 , 3 )  =-Binat  ( 5 ,  i )  *nc  ( i ) +omega  { 2 , 3 ) 

omega(l,3)=-Bmat(6,i) *nc(i)+omega(l,3) 

1  continue 
do  2  i=4, 6 

omega (1,2) =-Bmat ( 4 , i ) *nc { i ) +omega (1,2) 
omega (2,3) =-Bma t ( 5 , i ) *nc ( i ) +omega (2,3) 
omega(l,3)=-Bmat (6, i) *nc (i) +omega (1, 3 ) 

2  continue 

omega (2,1) =-omega (1,2) 
omega (3,2) =-omega (2,3) 
omega (3,1) =-omega (1,3) 
do  7  i=l,3 
do  8  j=l,3 

omega ( i , j ) =omega ( i , j ) /dsqr t ( 2 . OdO ) 
continue 
continue 

call  tenmatprod(Amat,nc, dine) 
c=1.000 
a=c/wil 
b=c/wi2 

if  (dabs (a-b) .gt. 0.01)  then 

omega (1,2) =omega (1 , 2 ) - (a*a+b*b) *dinc (4) / (dsqrt (2 . OdO) * (a*a’ 
endif 

if  (dabs(a-c)  .gt.0.01)  then  * 

omega(l,3)=omega(l,3) -(a*a+c*c) *dinc(6) / (dsqrt (2 . OdO) * (a*a- 
endif 

if  (dabs(c-b) .gt.0.01)  then 

omega (2,3) =omega (2,3)- (b*b+c  *c ) *dinc ( 5 ) / ( dsqrt ( 2 . OdO ) * (b*b 
endif 

omega (2,1) =-omega (1,2) 
omega (3,2) =-omega (2,3) 
omega (3,1) =-omega (1,3) 
omegac ( 1 ) =omega (1,1) 
omegac ( 2 ) =omega (2,2) 
omegac ( 3 ) =omega ( 3 , 3 ) 
omegac (4) =omega (1,2) *dsqrt (2 .DO) 
omegac ( 5 ) =omega (2,3) ^dsqr t ( 2 . DO ) 
omegac ( 6 ) =omega (1,3) *dsqr t ( 2 . DO ) 


s(l,l)=sc(l) 

s(2,2)=sc(2) 

s(3,3)=sc(3) 

s(l,2)=sc(4) /dsqrt (2. DO) 
s (2, 3) =sc(5) /dsqrt (2 .DO) 
s(3,l)=sc(6) /dsqrt (2. DO) 
s(2,l)=s(l,2) 
s(3,2)=s(2,3) 
s(l,3)=s(3,l) 

call  matprod(s, omega, dummy 1) 
call  matprod (omega,  s,durnmy2) 
do  300  i=l,3 
do  310  j=l,3 

sn  ( i ,  j  )  =duinmyl  ( i ,  j  )  -d\iinmy2  ( i ,  j  ) 

310  continue 
300  continue 

call  tenmatprod (Ce, nc, dummy) 
sn(l,l)=sn(l,l) -dummy (1) 
sn{2,2)=sn(2/2) -dummy (2) 
s n ( 3 , 3 ) = s n ( 3 , 3 ) - dummy ( 3 ) 
sn(l, 2) =sn(l, 2) -dummy (4) /dsqrt (2 .DO) 
sn(2, 3) =sn(2, 3) -dummy (5) /dsqrt (2 .DO) 
sn(3 , 1) =sn(3 , 1) -dummy (6) /dsqrt (2 .DO) 
sn(2,l)=sn(l,2) 
sn(3,2)=sn(2,3) 
sn (1 , 3) =sn (3 , 1) 
do  320  i=l,3 


do  330  j=l,3 
dummyl ( i , j ) =sn ( i , j ) 

330  continue 

320  continue 

c - STEP  3:  put  them  together 


snl (1,1) =sige(l, 1) +dlam*sn(l, 1) 
snl (2,2) =sige (2,2) +dlam*  sn ( 2 , 2 ) 
snl (3,3) =s ige (3,3) +dlam*  sn ( 3 , 3 ) 
snl(l,2)=sige(l,2)+dlam*sn(l,2) 
snl (2, 3) =sige(2, 3)+dlam*sn(2, 3) 
snl (3, 1) =sige(3, l)+dlam*sn(3, 1) 
snl (2,1) =snl (1,2) 
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snl(3,2)=snl(2,3) 
snl  (1, 3)  =snl  (3,1) 

These  components  of  snl  are  now  NOT  in  the  coordinate  farme 
coincides  with  the  orientation  of  the  inclusions.  This  is  due 
the  algorithm  used  (see  notes  for  more  details)  .  The  stress  coitponents 
are  now  relative  to  a  frame  which  is  obtained  by  rotating  the 
orientation  of  the  particles  by  R.  So  obtain  the  components  relative 
to  the  orientation  of  the  particles  at  the  beginning  of  this 
increment : 

call  matprod(R, snl,diainmyl) 
call  matprod ( dummy  1,RT,  snl) 


estimate  the  new  state  variables 
if  (evolf .eq. . true. )  then 
fn=f+(n(l,l)+n(2,2)+n(3,3) ) * (1.0-f) *dlam 
else 
fn=f 
endif 

if  (inpt.eq.l)  then 
endif 

call  Atensor  ( a, b,  c,  ad, bd,  cd,  f , mul ,  kl , mu2 ,  k2 ,  Amat) 
dummy  ( 1 )  =Amat  (3,1)  -Amat  (1,1) 
dummy  ( 2 )  =:Amat  (3,2)  -Amat  (1,2) 
dummy  ( 3 )  =Ainat  (3,3)  -Amat  (1,3) 
dummy ( 4 ) =Amat (3,4) -Amat (1,4) 
dummy ( 5 ) =Amat (3,5) -Amat (1,5) 
dummy ( 6 ) =Amat (3,6) -Amat (1,6) 
call  rowcolumnprod  (dummy,  dummy5 ,  d\im5 ) 
write (16, *) 'dummyl' , dummy 
write (16, *) 'dummy5' ,dummy5 
write ( 16 , * ) , ' A1 ' , Amat (3,1) 
write (16, *) , 'Al' , Amat (3,2) 
write (16, *) , ' A1 ' , Amat (3 , 3 ) 
write (16, *) , ' A1 ' , Amat (3 , 4) 
write (16, *) , ' A1 ^ Amat (3 , 5) 
write (16, *) , ' A1 ' , Amat (3,6) 
write (16,*) , ' A1 ' , Amat (1, 1) 
write (16,*) , ' A1 ' , Amat (1, 2 ) 
write (16, *) , ' A1 ' , Amat (1, 3 ) 
write ( 16 , * ) , ' A1 ' , Amat (1,4) 
write (16, *) , ' A1 ' , Amat (1 , 5) 
write (16,*) , 'Al' , Amat (1, 6) 
do  400  i=l,6 
dummyS ( i ) =nc ( i ) 
continue 

dummy ( 1 ) =Amat (3,1) -Amat (2,1) 
dummy ( 2 ) =Amat (3,2) -Amat (2,2) 
dummy ( 3 ) =Amat (3,3) -Amat (2,3) 
dummy  ( 4 )  =Ainat  (3,4)  -Amat  (2,4) 
dummy ( 5 ) =Amat (3,5) -Amat (2,5) 
dummy ( 6 ) -Amat (3,6) -Amat (2,6) 
call  rowcolumnprod  (dummy,  dummy 5,  dumlO) 
write  (16,  *)  'd\imitiy2  ' ,  dummy 
write  (16,  *)  'duminy52  '  ,d\immy5 
write (16, *) , ' A2 ' , Amat (3 , 1) 
write (16, *) , ' A2 ' , Amat (3,2) 
write (16, *) , ' A2 ' , Amat (3 , 3 ) 
write (16, *) , ' A2 Amat (3,4) 
write  (16,  *)  ,  '  A2  ' ,  Amat  (3, 5) 
write  (16,*),  '  A2  ' ,  Amat  (3 , 6) 
write  (16,  *) ,  'A2  ' ,  Amat  (2, 1) 
write(16,*),  'A2 ',  Amat  (2,2) 
write  (16,  *)  ,  'A2  ' ,  Amat  (2,3) 
write(16,*), ' A2 ' , Amat (2, 4) 
write(16,*), ' A2 ' , Amat (2 , 5) 
write (16, *) , ' A2 ' , Amat (2 , 6) 
if  (evolwi.eq. .true.)  then 
wiln=wil+dlam*wil*dum5 
wi  2n= wi  2 + dl  am*  wi  2  *  duml  0 
write  (15,  *)  'win'  ,wiln,wi2n 
write(15,  *)  'wil'  ,wil,wi2, dum5,duml0 
else 

wiln=wil 

wi2n=wi2 

endif 

if  ( (wiln.le.0.0) .or. (wi2n. le . 0 . 0) )  then 
write (15, *) 'Aspects  have  become  negative' 
call  flush (15) 
neg= . true . 
return 
endif 

“  estimating  new  orientations 
if  (evolti .eq. . true. )  then 
athetan=atheta- (omega (2, 3) *cp+omega (1,3) *sp) *dlam 


if  (dabs(wi2n~1.0) .lt.0.01)  then 
athetcui=0 . 0 
endif 

if  (st.ne.O.O)  then 

aphin=aphi-  ( (omega(2, 3)  *ss/st)  -omegad,  3)  *cs/st)  *dlam 
apsin=apsi+(-omega(l,2)+{ct/st) *{omega(2,3) *ss-omega (1, 3 ) *cs) ) * 
&  dlam 

else 

apsin=0.0 
aphin=0 . 0 
endif 

call  matprod(R, rot j , rotjn) 
else 

aphin=aphi 
apsin=apsi 
athetan=atheta 
do  27  i=l,3 
do  28  j=l,3 
rotjn(i, j ) =rotj (i,  j) 

28  continue 

27  continue 

endif 

ctn=dcos (athetan) 
s  tn=ds in ( athetan ) 
cpn=dcos (aphin) 
csn=dcos (apsin) 
spn=dsin (aphin) 
ssn=dsin ( apsin) 

rota (1, 1) =csn*cpn-ssn*ctn*spn  . 

rota (1,2) =-csn*spn-ssn*ctn*cpn 

rotad,  3)  =ssn*stn 

rota (2,1) =ssn*cpn+csn*ctn*spn 

rota (2,2) =-ssn*spn+csn*ctn*cpn 

rota (2,3) =-csn*stn 

rota (3,1) =stn*spn 

rota ( 3 , 2 ) =stn*cpn 

rota (3, 3) =ctn 

call  inatprod( rota,  rot jn,rotn) 
do  77  i=l,3 
do  88  j=l,3 
rotnt ( i , j ) =rotn ( j , i ) 


88  continue 

77  continue 

- - 

Q -  update  stress  and  state  variables: 


c  The  stress  must  now  be  expressed  w.r.t  the  new  orientation  of  the 

c  particles  so  that  ddsdde  can  be  calculated, 

c  First  we  take  it  back  w.r.t  the  GLOBAL  axes: 

call  matprod (rot, snl, dummy) 
call  matprod (dummy, rott, snl) 

c  and  then  express  it  in  w.r.t  the  new  orientation  of  the  particles, 

call  matprod(rotnt, snl, dummy) 
call  matprod (dummy, rotn, snl) 
sc(l)=snl(l,l) 
sc(2)=snl(2,2) 
sc(3)=snl(3,3) 
sc(4)=snl (1,2) *dsqrt(2.D0) 
sc(5)=snl (2,3) *dsqrt(2,D0) 
sc(6)=snl(3,l)*dsqrt(2.D0) 
f=fn 

wil=wiln 

wi2=wi2n 

c=1.000 

a=c/wi'l 

b=c/wi2 

ad=a 

bd=b 

cd=c 

*do  37  i=l,3 
do  38  j=l,3 
rotj (i, j)=rotjn(i, j) 

38  continue 

37  continue 

atheta=athetan 
aphi= aphin 


apsi=:apsin 
do  707  i=l,3 
do  808  j=l,3 
rot(i, j)=rotn(i, j) 
rott ( i , j ) =rotnt ( i ,  j  ) 
808  continue 

707  continue 

c -  updating  sig_yield 

do  501  i=l,3 


502 

501 


503 


C 

c 


392 

391 

c 

C-- 

C-- 

399 


401 


471 

470 


443 


410 


420 
c — 

C-- 


1002 

1001 


1004 

1003 


do  502  3=1,3 
deps { i , j ) =dlam*  n ( i , j ) 
continue 
continue 

duinl=:(deps(l,l)+deps(2,2)+deps(3,3) )  /3.0 
do  503  i=l,3 

deps  { i ,  i )  =deps  ( i ,  i )  -duitil 
continue 

dep_plas=deps ( 1 , 1 ) *  *2  +deps (2,2) *  *2+deps ( 3 , 3 ) *  *2 
dep_plas=dep_plas+2 .0* (deps (1, 2) **2+deps (2,3) **2+deps (3 , 1) 
dep_plas= (2 . 0/3 . 0) *dep_plas 
eps_plas  tic=eps_plast ic+dsqrt ( dep_plas ) 
sigyln=sigy0* ( (1 . 0+eps_plastic) ** (1 . 0/3 . 0) ) 
sigyl=sigyln 
sigyln=sigyl 


**2) 


call  Meff active  ( a , b,  c ,  ad,  bd,  cd,  f , mul ,  kll , itiu2 ,  k2 , MHS , path) 
do  391  i=l,6 
do  392  j=l,6 

MHS  (i ,  j  )  =3 . 0*inul*MHS  (i,  j  ) 
continue 
continue 

call  yieldtest (MHS , sc , sigyl , f , yc , y) 

CALCULATE  ddsdde 
calculate  N 
do  399  i=l,6 
duininy5  ( i )  =sc  ( i ) 
continue 

call  teninatprod(MHS,duittmy5,nc) 
do  401  i=l,6 
nc(i) =nc (i) / (1 .0-f ) 
nc(i)=2 .0*nc{i) /sigyl 
continue 
path=. false. 

call  Meffective{a,b,  c,ad,bd,  cd,  f  ,inule,  k:le,mu2e,  ek2e,Ce,patn) 
path= . true , 
do  470  i=l,6 
do  471  j=l,6 
test (i, j ) =Ce (i, j ) 
continue 
continue 

call  ludcmp ( test , 6, 6, indx, det) 
do  443  3=1,6 
det=det*test ( j , j ) 
continue 

if  (det. It. 0.0)  then 

write (15, *) 'Ce  is  wrong  during  plastic  step' 
call  flush(15) 
endif 

do  410  i=l,6 
duinmy5  ( i )  =nc  ( i ) 
continue 

call  tenmatprod(Ce, dummy 5, dummy) 
do  420  i=l,6 
dummyS ( i ) =nc ( i ) 
continue 

call  rowcolumnprod(duinny,dummy5,L) 

-  calculating  H 

'  step  A  calculating  dphi/df 
H=0.0 

path= . true , 

if  (evolf .eq. .true. )  then 

call  Mef fective(a,b, c,ad,bd,cd, f , mul, kll, mu2, k2, MHS, path) 
do  1001  i=l,6 
do  1002  j=l,6 
MHS ( i , 3 ) =3 . 0  *mul *MHS ( i , j ) 
continue 
continue 

call  yieldtest (MHS, sc, sigyl, f,yc, phi 1) 
f test=f+0 . 001*f 

call  Mef f ective ( a , b, c , ad , bd, cd, f tes t , mul , kll , mu2 ,  k2 , MHS , path) 
do  1003  i=l,6 
do  1004  j=l,6 
MHS (i, j ) =3 . 0*mul*MHS (i , 3 ) 
continue 
continue 

call  yieldtest (MHS, sc , sigyl , f test , yc , phi2 ) 
dphidf=(phi2-phil) / (ftest-f ) 

H=0.0 

H= -dphidf  *(1.0-f)*(nc(l) +nc ( 2 ) +nc ( 3 ) ) 
endif 

step  B  —  calculating  dphi/dwil 
if  (evolwi . eq . . true . )  then 


c 


path= . true . 

call  Meffective(a,b,c,ad,bd,  cd,  f  ,mul,  kll,mu2,k2,MHS,path) 

do  2001  i=l,6 
do  2002  j=l,6 
MHS  ( i ,  j  )  =3 . 0  *inul  *MHS  ( i ,  j  ) 
continue 
continue 

call  yieldtest (MHS,sc, sigyl, f ,yc,phil) 

wiltest=wil+0 . 001*wil 

ate=c/wiltest 

call  Mef f ective (ate,b, c, ate, b,  c,  f  ,inul, kll ,inu2 ,  k2 , 

^  MHS, path) 

do  2003  i=l,6 
do  2004  j=l,6 
MHS  (i,  j  )  =3 . 0*inul*MHS  (i,  j  ) 
continue 
continue 

call  yieldtest (MHS,sc,sigyl, f ,yc,phi2) 
dphidwil= (phi2-phil) / (wiltest-wil) 

H=H-dphidwil  *wil  *duiti5 
endif 

-  step  C  —  calculating  dphi/dwi2 
if  (evolwi .eq. . true. )  then 
path=.true. 

call  Mef  f ective  (a,b,  c,ad,bd,  cd,  f,mul ,  kll, mu2 ,  k2, MHS, path) 
do  3001  i=l,6 
do  3002  j=l,6 
MHS ( i , j ) =3 . 0  *mul *MHS ( i , j ) 
continue 
continue 

call  yieldtest (MHS, sc, sigyl, f , yc,phil) 

wi2test=wi2+0 . 001*wi2 

bte=c/wi2test 

call  Meffective{a,bte,c,a,bte,c, f ,mul,kll,mu2,k2, 
i  MHS, path) 

do  3003  i=l,6 
do  3004  j=l,6 
MHS  ( i ,  j )  =3 . 0*inul*MHS  (i ,  j  ) 
continue 
continue 

call  yieldtest (MHS, sc, sigyl , f,yc,phi2) 
dphidwi2=(phi2-phil) / (wi2test-wi2 ) 

H=H-dphidwi2  *wi2  *duml0 
endif 

if (debug . eq. . true. )  then 

write { 15 , * ) ' updates ' 

call  flush (15) 

endif 

hrd=H 


L=L+H 

write(15,*) 'L',L,  'H',H 
call  flush(15) 
if  (L. It. 0.0)  then 
write (15, *) 'L  is  negative' 
call  flush(15) 
endif 

s(l,l)=sc(l) 

s(2,2)=sc(2) 

s(3,3)=sc{3) 

s(l#2)=sc(4) /dsqrt(2.D0) 
s(2,3)=sc(5) /dsqrt(2.D0) 
s(3,l)“sc{6) /dsqrt(2.D0) 
s(2,l)=s{l,2) 
s{3,2)=s(2,3) 
s(l,3)=s(3,l) 
do  450  i=l,6 
duinmy5  ( i )  =nc  ( i ) 
continue 

call  tenitiatprod(Ce, dummy 5 , dummy) 
call  columnrowprod( dummy, dummy, ddsdde) 
do  430  i=l,6 
do  440  j=:l,6 

ddsdde ( i, j ) =Ce (i, j ) -ddsdde (i, j ) /L 
continue 
continue 


do  111  i=l,6 
do  112  j=l,6 
test ( i , j ) =ddsdde ( i , j ) 
continue 
continue 

call  ludcmp ( test , 6 , 6 , indx, det) 


113 

c 


ccc 

422 

421 


4432 

4431 

c - 


571 

570 


543 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


do  113  j=l,6 
det=det*test (j, j) 
continue 

if  (det.lt. 0.0)  then 
write(15,*) 'detl',det 
call  flush (15) 
endif 

call  inatpr9d(S/ omega, duntmyl) 
call  matprod (omega,  s,duinmy2) 

do  421  i=l,3 
do  422  j=l,3 

dxoinmyl  ( i ,  j  )  =duinmyl  ( i ,  j  )  -duinmy2  ( i ,  j ) 
dummyl (i, j ) =0 . 0 
continue 
continue 

dummy  5  ( 1 )  =duinmyl  (1,1) 

dummy 5  ( 2 )  =dummyl  (2,2) 

dummy 5 ( 3 ) =dummyl (3,3) 

dummy5 ( 4 ) =dummyl (1,2) *dsqr t ( 2 . DO ) 

dummy5 ( 5 ) =dummyl (2,3) *dsqrt (2 . DO ) 

dummy 5 ( 6 ) =dummyl (3,1) *dsqrt ( 2 . DO ) 

call  CO lumnrowprod( dummy 5, dummy, ddsdum) 
if  (evolf .eq. . true. )  then 
do  4431  i=l,6 
do  4432  j=l,6 

ddsdde ( i , j ) =ddsdde ( i , j ) +ddsdum ( i , j ) /L 
continue 
continue 
endif 

Express  components  of  ddsdde  w.r.t  fixed  GLOBAL  axes, 
call  mat 2 tensor (ddsdde, dumm) 
call  rot 4order  (dumm,  rotn,  dumml) 
call  ten2matrix (dumml, ddsdde) 
do  570  i=l,6 
do  571  j=l,6 
dumd ( i , j ) =ddsdde ( i , j ) 
continue 
continue 

call  ludcmp (dumd, 6, 6, indx,det) 
do  543  j=l,6 
det=det*dumd( j  ,  j  ) 
continue 

if  (det.lt. 0.0)  then 
write (15 , *) 'abc' , a,b, c 
write ( 15 , * ) ddsdde 
write (15 , *) '  ' 

write(15, *) 'Ce' 
write(15, *)Ce 
write(15, *) 'det' ,det 
endif 

Express  stress  components  w.r.t  GLOBAL  axes, 
if  (evolti.eq. .true. )  then 
call  matprod(rotn, s, dummyl) 
call  ma tprod( dummyl, rotnt, s) 
endif 

sc(l)=s(l,l) 

sc(2)=s(2,2) 

sc(3)=s(3,3) 

sc(4)=s(l,2) *dsqrt(2.D0) 
sc(5)=s(2,3)*dsqrt(2.D0) 
sc(6)=s(3,l) *dsqrt(2.D0) 

If  the  constitutive  equation  is  written  in  terms  of  the 
Jaiimann  rate  of  the  Cauchy  stress,  then  ddsdde  will  be 
unsymmetric.  The  additional  quantity  that  appears  in  ddsdde 
is  added  below. 

write (15, *) 'ddsll ' , ddsdde (1,1) ,s(l,l) 
write (15, *) 'dds22 ' , ddsdde (2,2) ,s(2,2) 
write(15,*) 'dds33 ' , ddsdde (3 , 3) , s (3, 3) 
write (15, *) 'dds44 ', ddsdde (4, 1) ,sc(4) 
write (15, *) 'dds55' , ddsdde (5,1) ,sc(5) 
write(15, *) 'dds66' , ddsdde (6, 1) ,sc(6) 
call  flush (15) 

ddsdde(l,l)=ddsdde(l,l)+s (1,1) 
ddsdde(l,2)=ddsdde(l,2)+s(l,l) 
ddsdde(l,3)=ddsdde(l,3)+s(l,l) 
ddsdde (2,1) =ddsdde (2,l)+s(2,2) 
ddsdde (2,2) =ddsdde (2 , 2 ) +s (2 , 2 ) 
ddsdde (2,3) =ddsdde ( 2 , 3 ) +s (2 , 2 ) 
ddsdde (3,1) =ddsdde (3,l)+s(3,3) 
ddsdde (3 , 2) =ddsdde (3 , 2 ) +s ( 3 , 3 ) 
ddsdde (3,3) =ddsdde ( 3 , 3 ) +s { 3 , 3 ) 
ddsdde (4,1) =ddsdde (4,1) +sc ( 4 ) 
ddsdde(4,2)=ddsdde(4,2)+sc(4) 
ddsdde(4, 3)=ddsdde(4, 3)+sc(4) 


c 

c 

c 

c 

c 

c 

c~— 
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478 

477 


474 
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475 


531 

530 


538 
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c 


542 

541 


999 


2 
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ddsdde (5,1) =ddsdde (5,1) +sc ( 5 ) 
ddsdde(5,2)=ddsdde(5,2)+sc(5) 
ddsdde(5,3)=ddsdde(5,3)+sc(5) 
ddsdde(6,l)=ddsdde(6,l)+sc(6) 
ddsdde(6, 2)=ddsdde(6,2)+sc(6) 
ddsdde ( 6 , 3 ) =ddsdde ( 6 , 3 ) +SC ( 6 ) 

Interchange  con^onents  of  ddsdde  here  (since  our  notation  for 
stress  and  strain  are  different  from  what  they  use) . 
do  477  i=i;3 
ao  478  j=4,6 

ddsdde ( i , j ) =ddsdde ( i , j ) /dsqrt ( 2 . DO ) 
continue 
continue 
do  473  i=4,6 
do  474  j=l,3 

ddsdde ( i , j ) =ddsdde ( i , j ) /dsqrt ( 2 . DO ) 
continue 
continue 
do  475  i=4,6 
do  476  j=4,6 

ddsdde ( i , j ) =ddsdde ( i , j ) / ( 2 . DO ) 
continue 
continue 
do  530  i=l,6 
do  531  j=l,6 
dumd ( i , j ) =ddsdde ( i , j ) 
continue 
continue 
do  538  i=l,4 
ddsdde ( i , 5 ) =dumd ( i , 6 ) 
ddsdde ( i , 6 ) =dumd ( i , 5 ) 
continue 
do  539  i=l,4 
ddsdde ( 5 , i ) =dumd ( 6 , i ) 
ddsdde ( 6 , i ) =dumd ( 5 , i ) 
continue 

ddsdde (5,5) =dumd (6,6) 
ddsdde (5,6) =dumd (6,5) 
ddsdde (6,5) =dumd (5,6) 
ddsdde (6,6) =dumd (5,5) 

-  The  above  is  the  final  ddsdde 
do  541  i=l,6 
do  542  j=l,6 
diimd  ( i ,  j  )  =ddsdde  ( i ,  j  ) 
continue 
continue 

if (debug. eq. . true. )  then 

write ( 15, * ) 'updatelO ' 

call  flushdS) 

endif 

return 

end 


subroutine  stre(ss, sse,p,q, estimate, Ce, nine) 
real*8  MHS(6,6) 

real*8  p(6),q(6),r(6) , eps_plastic,Ce (6, 6) 

real* 8  mul,kl,mu2,k2, sigyl, f ,wil,wi2, sc (6) ,strainc(6) 

real *8  a, b, c, ad,bd, cd, sigyO, kll ,y,hrd, ss (6) , sse (6) , estimate (6) 

integer  nine , i , j , ic 

logical  path, check, yield, neg, yc , debug 
logical  evolf , evolwi , evolti 
common  /controls/evolf, evolwi, evolti 
common*  /mkmodulus/mul, kl,mu2, k2 , sigyO , kll 

common  /mkdatal / f , wil , wi2 , sc , strainc , eps^plas tic , sigyl , hrd 

common  /mkdata2 /yield, check, neg 

common  /mdebug/ debug 

path= . true , 

c=1.000 

a=c/wil 

b=c/wi2 

ad=a 

bd=b 

cd=c 

call  Mef f ective (a,b, c, ad, bd, cd, f ,mul , kll,mu2 , k2 ,MHS, path) 
do  1  i=l,6 
do  2  j=l,6 

MRS (i, j ) =3 . 0*mul*MHS (i , j ) 
continue 
continue 

i f  ( check . eq . . true . )  then 
call  yieldtest (MHS, ss, sigyl, f,yc,y) 
if  (yc.eq. . true. )  then 
write (15, *) 'WARNINGl  IN  STRE',y 


stop 
endif 

call  yieldtest  (MHS,  sse,  sigyl,  f , yc,y) 
if  (yc.eq, .false. )  then 
write (15, *) 'WARNING2  IN  STRE' 
stop 
endif 
endif 

do  100  ic=i,ninc 
do  10  i=l,6 

r{i)=0.5*(p(i)+q(i)) 

10  continue 
call  tenmatprod{Ce,r, estimate) 

do  11  i=l,6 

estimate ( i ) =es timate ( i ) +ss ( i ) 

11  continue 

call  yieldtest (MHS, estimate, sigyl, f,yc,y) 
cwrite(15,*) 'est  in  stre',y 

if  ( (yc.eq. .true. ) .and. (y.lt .1 . Od-10) ) then 
return 
else 

if  (yc. eq. . false . )  then 
do  20  i=l,6 
p(i)=r(i) 

20  continue 

else 

do  30  i=l,6 
q(i)=r(i) 

30  continue 

endif 
endif 

100  continue 

writedS,  *) 'problem  in  stre' 
stop 
return 
end 

Q -  This  is  used  to  decompose  the  F  tensor, 

subroutine  decompose (F,R,U, C,D) 
real*8  F (3, 3) ,R(3, 3) ,FT(3, 3) ,C(3 , 3) 
real*8  D(3) ,V(3,3) ,U(3,3) 
integer  n,np, nrot, i, j 
logical  debug 
common  /mdebug/ debug 
do  10  i=l,3 
do  20  j=l,3 
FT(i, j)=F(j,i) 

20  continue 

10  continue 

call  matprod(FT,F,C) 
n=3 
np=3 

call  jacobi {C,n,np,D,V,nrot) 
do  30  i=l,3 
do  40  j=l,3 
C(i,j)=0.0 
U(i, j)=0.0 
40  continue 

30  continue 

do  50  i=l,3 
do  60  j=l,3 
FT(i,  j)=V(i,l)*V(j,l) 

U  ( i-,  j  )  =FT  ( i ,  j  )  *dlog  ( dsqrt  (d  ( 1 )  ) ) 

C(i, j)=FT(i, j) /dsqrt (d(l) ) 

60  continue 

50  continue 

do  70  i=l,3 
do  80  j=l,3 
FT(i,  j)=V(i,2)*V(j,2) 

U(i,j)  =  (FT(i,j)  *dlog  (dsqrt  (d  (2 ) )  ) )  +U  ( i ,  j ) 

C(i, j)=(FT(i, j)/dsqrt(d(2) ) )+C(i, j) 

80  continue 

70  continue 

do  90  i=l,3 
do  100  j=l,3 
FT(i,  j)=V(i,3)*V(j,3) 

U(i,j)=(FT(i,j) *dlog (dsqrt (d(3) ) ) ) +U(i, j) 

C(i, j)=(FT(i,j) /dsqrt (d(3) ) )+C(i, j) 
continue 
continue 

call  matprod (F, C,R) 
call  matprod (R,V, C) 

C(1,1)=C(1,1) /dsqrt (C{ 1,1) **2+C(2,l) **2+C ( 3 , 1 ) **2 ) 


100 

90 


C(2,1)=C(2,1) 

C(3,1)=C(3,1) 

C(1,2)=C(1,2) 

C(2,2)=C(2,2) 

C(3,2)=C(3,2) 

C(1,3)=C(1,3) 

C(2,3)=C(2,3) 

C(3,3)=C(3,3) 

return 

end 


/dsqrt(C(l, 
/dsqrt(C(l, 
/dsqrt(C(l, 
/dsqrt (C(l, 
/dsqrt(C(l, 
/dsqrt (C(l, 
/dsqrt (C(l, 
/dsqrt (C(l, 


1)**2+C(2,1) 

1) **2+C(2,1) 

2) **2+C(2,2) 
2)**2+C{2,2) 

2)  **2+C(2,2) 

3) **2+C(2,3) 
3) **2+C(2,3) 
3)**2+C(2,3) 


**2+C(3,l) **2) 
**2+C(3,l) **2) 
**2+C(3,2) **2) 
**2+C(3,2) **2) 
**2+C(3,2) **2) 
**2+C(3,3) **2) 
**2+C(3,3)**2) 
**2+C(3,3)**2) 


c _  This  is  used  to  decompose  the  F  tensor . 

subroutine  decomposel (F, R,U) 
real*8  F (3, 3) ,R{3, 3) ,FT(3 , 3) ,C(3 , 3) 
real*8  D(3) ,V(3,3) ,U(3,3) ,Uinv{3,3) 
integer  n,np,nrot, i, j 
logical  debug 
common  /mdebug /debug 
do  10  i=l,3 
do  20  j=l,3 
FT(i, j)=F(j,i) 

20  continue 

10  continue 

call  matprod(FT,F,C) 

n=3 

np=3 

call  jacobi (C,n,np,D,V,nrot) 
do  30  i=l,3 
do  40  j=l,3 
C(i, j)=0.0 
40  continue 

30  continue 

do  50  i=l,3 
do  60  j=l,3 
FT(i, j)=V(i,l)*V(j,l) 

U(i,  j ) =FT (i, j ) * (dsqrt (d(l) ) ) 

C(i,  j )=FT(i, j) /dsqrt (d(l) ) 

60  continue 

50  continue 

do  70  i=l,3 
do  80  j=l,3 
FT(i, j)=V(i,2)*V(j,2) 

U(i, j)=(FT(i, j)*(dsqrt(d(2) ) ) )+U(i, j) 
C(i,  j)  =  (FT(i, j)/dsqrt(d(2) ) )+C(i,j) 

80  continue 

70  continue 

do  90  i=l,3 
do  100  j=l,3 
FT{i, j)=V{i,3)*V(j,3) 

U(i, j)=(FT(i, j)*(dsqrt(d(3) ) ) )+U(i, j) 
C(i, j)={FT{i, j)/dsqrt(d(3) ) )+C(i, j) 
100  continue 

90  continue 

do  110  i=l,3 
do  120  j=l,3 
V(i, j)=U(i, j) 

120  continue 

110  continue 

call  inverse3x3 (V,Uinv) 
call  inatprod(F,Uinv, R) 
return 
end 


c 


C _ This  is  a  routine  to  multiply  2nd  order  tensors  (3x3  matrices) 

subroutine  matprod (p, q, r) 
real*8  p(3,3) ,q(3, 3) , r (3 , 3) 
integer  i,j,k 
‘  do  10  i=l,3 
do  11  j=l,3 
r(i,j)=0.0 

11  continue 

10  continue 

do  12  i=l,3 
do  13  j=l,3 
do  14  k=l,3 

r(i, j)=r(i, j)+p(i,k) *q(k, j) 

14  continue 

13  continue 

12  continue 
return 
end 


SUBROUTINE  JACOBI (A, N, OT, D,V,NROT) 
implicit  double  precision  (a-h,  o-z) 
integer  niriax 
PARAMETER  (NMAX=100) 

real*8  A(NP,NP) ,D(NP) ,V(NP,NP) ,B(NMAX) , Z (NMAX) 
integer  ip, iq, i,n, np.nrot, j 
DO  12  IP=1,N 
DO  11  IQ=1,N 
V(IP,IQ)=0. 

11  CONTINUE 
V(IP,IP)=1. 

12  CONTINUE 

DO  13  IP=1,N 
B(IP)=A(IP, IP) 

D(IP)=B(IP) 

Z(IP)=0. 

13  CONTINUE 
NROT=0 

DO  24  1=1,50 
SM=0. 

DO  15  IP=1,N-1 
DO  14  IQ=IP+1,N 

SM=SM+DABS(A(IP, IQ) ) 

14  CONTINUE 

15  CONTINUE 
IF(SM.EQ.O. ) RETURN 
IF{I.LT.4)THEN 

TRESH=0 . 2  *SM/N**2 
ELSE 

TRESH=0 . 

ENDIF 

DO  22  IP=1,N-1 
DO  21  IQ=IP+1,N 

G=100.*DABS(A(IP,IQ) ) 

IF( (I.GT.4) .AND. (DABS(D(IP) ) +G.EQ.DABS (D (IP)  )  ) 
*  .AND. (DABS(D(IQ) ) +G.EQ.DABS (D(IQ) ) ) ) THEN 

A(IP,IQ)=0. 

ELSE  IF(DABS(A(IP,IQ) ) .GT.TRESH)THEN 
H=D(IQ)-D(IP) 

IF (DABS (H)+G.EQ. DABS (H) ) THEN 
T=A(IP,IQ)/H 
ELSE 

THETA=0 . 5*H/A(IP, IQ) 

T=1 . / (DABS (THETA) +DSQRT ( 1 . DO+THETA*  *2 ) ) 

IF (THETA.lt. 0 . ) T=-T 
ENDIF 

C=1 . /DSQRT ( 1 . DO+T*  *2 ) 

S=T*C 

TAU=S/ (l.+C) 

H=T*A(IP, IQ) 

Z(IP)=Z(IP)-H 

Z(IQ)=Z(IQ)+H 

D(IP)=D(IP)-'H 

D(IQ)=D(IQ)+H 

A(IP,IQ)=0. 

DO  16  J=1,IP-1 
G=A(J,IP) 

H=A(J,IQ) 

A(  J, IP) =G-S* (H+G*TAU) 

A( J, IQ) =H+S* (G-H*TAU) 

16  CONTINUE 

DO  17  J=IP+1,IQ-1 
G=A(IP, J) 

-  H=A(J,IQ) 

A(IP,J)=G-S* (H+G*TAU) 

A(J, IQ) =H+S* (G-H*TAU) 

17  CONTINUE 

DO  18  J=IQ+1,N 
G=A(IP, J) 

H=A(IQ,J) 

A(IP,J)=G-S*(H+G*TAU) 

A(IQ, J) =H+S* (G-H*TAU) 

18  CONTINUE 

DO  19  J=1,N 
G=V(J,IP) 

H=V(J,IQ) 

V(J, IP) =G-S* {H+G*TAU) 

V(J,IQ)=H+S* (G-H*TAU) 

19  CONTINUE 
NROT=NROT+l 

ENDIF 

21  CONTINUE 

22  CONTINUE 

DO  23  IP=1,N 


o  O  O  CN  rH 

CNrHU  U  rororoU 


B(IP)=B{IP)+Z(IP) 

D(IP)=:B{IP) 

Z(IP)=0. 

23  CONTINUE 

24  CONTINUE 

PAUSE  '50  iterations  should  never  happen' 

RETUKN 

END 


C - This  is  to  calculate  the  A  tensor 

subroutine  Atensor (a,b/ c, ad,bd, cd, c2 ,mul, kl,niu2 , k2 , Amat) 
real *8  Amat (6, 6) ,Li (6, 6) ,L2 (6/ 6) , Si (6, 6) 
real*8  Ml (6, 6) ,mul,mu2,kl,k2,nul,nu2 
real*8  a,b, c, ad,bd, cd, cl, c2 

real*8  d(6,6) , dl (6, 6) , d2 (6, 6) ,pil212 ,pi2323 ,pil313 

integer  i,j 

logical  debug, spath 

common  /paths /spath 

common  /mdebug  / debug 

spath=.true. 

nul=0.5* (3 .0*kl-2 . 0*mul) / (3 . 0*kl+mul) 
nu2=0. 5*  (3 . 0*k2*-2 . 0*mu2)  /  (3 .0*k2+mu2) 
cl=1.0*-c2 

if (debug. eq. . true. )  then 
write (15, *) ' Atenl ' 
call  flush(15) 
endif 

c  Initialize  the  modulus  tensors 

do  10  i=l,6 
do  20  j=l,6 
LKi,  j)=0.0 
L2(i, j)=0.0 
d(i, j)=0.0 
continue 
continue 

Obtain  the  (6x6)  modulus  tensors,  LI  and  L2 
Ll(l,l)=kl+(4.0/3.0) *mul 
Ll(2,2)=kl+(4.0/3.0) *mul 
Ll(3,3)=kl+(4.0/3.0) *mul 
Ll(l,2)=:kl-(2.0/3.0)  *mul 
LI (2,l)=kl-(2.0/3.0) *mul 
Ll(2,3)=kl-(2.0/3.0) *mul 
Ll(3,2)=kl-(2.0/3.0) *mul 
Ll(3,l)=kl- (2. 0/3.0) *mul 
LI ( 1 , 3 ) =kl- (2 . 0/3 . 0)  *mul 
Ll(4,4)=2.0*mul 
Ll(5,5)=2.0*mul 
Ll(6,6)=2.0*mul 
L2(l,l)=k2+(4.0/3.0) *mu2 
L2(2,2)=k2+(4.0/3.0) *mu2 
L2(3,3)=k2+(4.0/3.0)*mu2 
L2(l,2)=k2-(2.0/3.0) *mu2 
L2(2,l)=k2-(2.0/3.0)*mu2 
L2 (2,3)=k2-(2.0/3.0)*mu2 
L2(3,2)=k2-(2.0/3.0) *mu2 
L2(3,l)=k2-(2.0/3.0)*mu2 
L2(l,3)=k2-(2.0/3.0) *mu2 
L2(4,4)=2.0*mu2 
L2(5,5)=2.0*mu2 
L2(6,6)=2.0*mu2 
Begin  calculating  A  tensor 
if (debug. eq. . true. )  then 
write  (1-5,  *)  '  Aten2  ' 
call  flushdS) 
endif 

call  inverse (LI, Ml) 
call  tenprod(Ml,L2,d) 
do  30  i=l,6 
d(i,i)=d(i,i)-1.0 
continue 
do  31  i=l,6 
do  32  j=l,6 
dl(i, j)=0.0 
d2(i, j)=0.0 
continue 
continue 

calculate  S  tensors 
if (debug. eq. . true. )  then 
write ( 15 , * ) ' Aten3 ' 
call  flush(15) 
endif 

call  stensor (a,b, c,nul ,kl,mul, Si,pil212,pil313 ,pi2323) 
if (debug. eq. . true. )  then 


write { 15 , * ) ' Aten4 ' 
call  flush (15) 
endif 

do  60  i=l,6 
do  70  j=l,6 

dl(i, j)=Si(i, j)-c2*Si(i, j) 
70  continue 
60  continue 

call  tenprod(dl, d, d2) 
do  80  i=l,6 
d2(i,i)=d2(i,i)+1.0 
80  continue 

if (debug . eq . . true . )  then 
write (15,*) 'Aten5',d2 
call  flush (15) 
endif 

call  inverse  (d2  ,Aitiat) 

if (debug. eq. .true. )  then 

write(15, *) 'Aten6' ,  Amat 

call  flush (15) 

endif 

return 

end 


subroutine  columnrowprod (a,b, c) 
real*8  a(6) ,b(6) ,c (6, 6) 
integer  i,j 
do  10  i=l,6 
do  20  j=l,6 
c(i, j)=a(i) *b(j) 

20  continue 
10  continue 
return 
end 


C - This  is  to  calculate  the  effective  modulus  of  the  composite. 

subroutine  Mef fective {a,b, c,ad,bd, cd, c2,mul,kl,mu2,k2,MHS,path) 

real*8  LHS{6, 6) ,L1 (6, 6) ,L2 (6, 6) , Si (6, 6) ,MHS(6, 6) 

real *8  Ml (6, 6) ,mul,mu2 ,kl, k2,nul,nu2 

real*8  a, b, c, ad,bd, cd, cl, c2 

real*8  d(6, 6) ,dl (6, 6) ,pil212,pi2323,pil313 

integer  i,j 

logical  path, spa th, debug 
common  /paths/spath 
common  /mdebug/ debug 
spath= . true. 

if (debug. eq. . true . )  then 
write(15, *) 'Meffl' 
call  flush(15) 
endif 

nul=0 . 5*  (3 .0*kl-2 . 0*mul)  /  (3 .0*kl+mul) 
nu2=0.5* (3.0*k2-2.0*mu2) / (3.0*k2+mu2) 
cl=1.0-c2 

c  Initialize  the  modulus  tensors 

do  10  i=l,6 
do  20  j=l,6 
Ll(i, j)=0.0 
L2(i, j)=0.0 
20  continue 
10  continue 

c  Obtain  the  (6x6)  modulus  tensors,  LI  and  L2 

Ll(l,l)=kl+(4.0/3.0) *mul 
Ll(2,2)=kl+(4.0/3.0) *mul 
Ll(3,3)=kl+(4.0/3.0)*mul 
Ll(l,2)=kl-(2.0/3.0)*mul 
Ll(2,l)=kl-(2.0/3.0) *mul 
Ll(2,3)=kl-(2.0/3.0)*mul 
Ll(3.2)=kl-(2.0/3.0) *mul 
Ll(3,l)=kl-(2.0/3.0)*mul 
Ll{l,3)=kl-(2.0/3.0)*mul 
Ll(4,4)=2.0*mul 
Ll(5,5)=2,0*mul 
Ll(6,6)=2.0*mul 
L2 (l,l)=k2+(4.0/3.0) *mu2 
L2 (2,2)=k2+(4.0/3.0) *mu2 
L2 (3 , 3 ) =k2+ (4 . 0/3 . 0) *mu2 
L2(l,2)=k2-(2.0/3.0)*mu2 
L2(2,l)=k2-(2.0/3.0)*mu2 
L2 (2,3)=k2-(2.0/3.0) *mu2 
L2 (3,2)=k2-(2.0/3.0) *mu2 
L2 (3,l)=k2-(2.0/3.0) *mu2 
L2 (1, 3)=k2-(2 .0/3 .0) *mu2 


L2(4,4)=2.0*mu2 
L2(5,5)=2.0*mu2 
L2{6,6)=2.0*inu2 
c  Begin  calculating  LHS 

do  21  1=1,6 
do  22  j=l,6 
dl(i, j)=Ll(i, j) 

22  continue 
21  continue 

i f ( debug . eq . . true . )  then 
writedS,*)  'Meff2' 
call  flushdS) 
endif 

call  inverse (dl, Ml) 
call  tenprod(Ml,L2,d) 
i f ( debug . eq . . true . )  then 
write (15 , *) 'Mef f 3 ' 
call  flushdS) 
endif 

do  30  i=l,6 
d(i,i)=d(i,i)-1.0 

30  continue 

do  31  i=l,6 
do  32  j=l,6 
did,  j)=0.0 
32  continue 

31  continue 

call  inverse (d,dl) 
if (debug . eq. . true. )  then 
writedS,  *)  'Meff4' 
call  flushdS) 
endif 

c  calculate  S  tensors 

if (debug. eq. . true . )  then 
write (15, *) 'MeffS' ,a,b, c 
call  flushdS) 
endif 

call  s tensor (a, b, c , nul , kl,mul , Si,pil212 ,pil313 ,pi2323) 
if (debug. eq. . true , )  then 
writedS,*)  'Meff6' 
call  flushdS) 
endif 

if (debug . eq. . true . )  then 
writedS,  *)  'Meff6a' 
call  flushdS) 
endif 

do  40  i=l,6 
do  50  j=l,6 
d(i, j)=0.0 
50  continue 
40  continue 

i f ( debug . eq . . true . )  then 
writedS,  *)  'Mef  f  6a_l ' ,  Si 
call  flushdS) 
endif 

do  60  i=l,6 
do  70  j=l,6 

if (debug. eq. . true . )  then 

writedS,  *)' j  =  ' ,  j  /  '  i=' ,  i/ Si  (i,  j  ) 

call  flushdS) 
endif 

d(i, j ) =Si (i, j ) -c2*Si (i, j ) 

70  continue 

60  continue 

if (debug. eq. . true . )  then 
writedS,  *)  'Meff6b' 
call  flushdS) 
endif 

do  80  i=l,6 
do  90  j=l,6 
d(i, j)=dl(i, j)+d(i, j) 

90  continue 

80  continue 

if (debug. eq. . true. )  then 
write(15,*) 'Meff6c' 
call  flushdS) 
endif 

do  100  i=l,6 
do  110  j=l,6 
did,  j)=0.0 
110  continue 

100  continue 

if (debug. eq. . true . )  then 
write (15 , * ) 'Mef f 7 ' 


call  flush (15) 
endif 

call  inverse (d,dl) 
if (debug. eq. . true. )  then 
writedS,*)  'MeffTa' 
call  flush (15) 
endif 

do  120  i=l,6 
do  130  j=l,6 
d(i, j)=c2*dl(i,  j) 

130  continue 
120  continue 

do  140  i=l,6 
d(i,i)=d(i,i)+1.0 
140  continue 

do  150  i=l,6 
do  160  j=1.6 
dl(i, j)=0.0 

160  continue 
150  continue 

call  tenprod(Ll, d, LHS) 
do  161  i=l,6 
do  162  j=l,6 
d(i, j)=0.0 
d(i, j)=LHS(i, j) 

162  continue 

161  continue 

call  inverse (d,MHS) 
if  (path. eq. . false . )  then 
do  170  i=l,6 
do  180  j=l,6 
MHS(i,3)=LHS(i, j) 

180  continue 

170  continue 

endif 

c - symmetrize  M  to  avoid  small  numerical  variations 

do  190  i=l,6 
do  200  j=l,6 
Ml(i, j)=MHS(i, j) 

200  continue 

190  continue 

do  210  i=l,6 
do  220  j=l,6 

MHS(i, j)=C-5*(Ml (i, j)+Ml(j,i)  ) 

220  continue 

210  continue 

if (debug. eq. . true. )  then 
write(15,*) 'Meff9' 
call  flush(15) 
endif 
return 
end 


C _  This  is  to  evaluate  the  S  tensor  for  an  ellipsoidal  inclusion. 

subroutine  s tensor (aold, bold, cold, nu , kl , mul , s , pi 2 12 , pl3 13 , p2 323 ) 
real*8  ia, ib, ic, iaa, ibb, icc, iab, ibc, ica, iac, icb, iba 
real*8  a,b,c,p,el2 
real *8  ff,e 

real*8  sllll , s2222 , s3333 , sll22 , sll33 , s2233 , s2211, s3311 , s3322 

real*8  sl212 , sl313 , s2323 

real*8  q, r , nu, kl ,mul 

real *8  x, qqc, qqc2, aa,bb 

real*8  s(6,6) 

real*8  al,bl,cl, si (6, 6) , aold, bold, cold 
real*8  pl212,pl313,p2323,ppl212,ppl313,pp2323 
real*8  pl221,pl331,p2332,ppl221,ppl331,pp2332 
real*8  p2112,p3113,p3223,pp2112,pp3113,pp3223 
real*8  p2121,p3131,p3232,pp2121,pp3131,pp3232 
integer  ip,i,j 
logical  spath, debug 
common  /paths /spath 
common  /mdebug / debug 
p=3. 14159265358979 
pl212=0.0 
pl313=:0.0 
p2323=0.0 
p=3. 14159265 

if  (dabs (aold-1. DO) .It. 0.001)  then 
a=l.D0 
else 
a=aold 
endif 


c 


if  (dabs (bold-1. DO) .It. 0.001)  then 
b=l.DO 
else 
b=bold 
endif 

if  (dabs(cold-l.DO) .lt.0.001)  then 
c=l.D0 
else 
c=cold 
endif 
do  i=l, 6  . 
do  j=l,6 
s(i, j)=0.0 
sl(i, j)=0.0 
enddo 
enddo 
ip=0 
ia=0 . 0 
ib=0 . 0 
ic=0 . 0 
iab=0 . 0 
ibc=0 . 0 
iac=0 . 0 
iaa=0.0 
ibb=0.0 
icc=0.0 

if  ( (a.eq.b) .and. (b.eq.c) )  then 

sllll=(8.D0*miil+9.D0*kl) /(5.D0* (4 .D0*mul+3 .D0*kl) ) 
s2222=(8 .D0*mul+9 .D0*kl) / (5 .DO*  U .D0*mul+3 . D0*kl) ) 
s3333=(8.D0*mul+9.D0*kl) /(5.D0* (4 .D0*mul+3 .D0*kl) ) 
sll22=(3 .D0*kl-4.D0*mul) /(5.D0* (4 .D0*mul+3 .D0*kl) ) 
sll33=(3.D0*kl-4.D0*mul) /(5.D0* (4 .D0*mul+3 .D0*kl) ) 
s2233=(3 .D0*kl-4.D0*mul) / (5. DO* (4 .D0*mul+3 .D0*kl) ) 
s2211=(3.D0*kl-4.D0*mul) /(5.D0* (4 .D0*mul+3 .D0*kl) ) 
s3311=(3 .D0*kl-4.D0*mul) / (5 .DO* (4 .D0*mul+3 .D0*kl) ) 
s3322= (3 .D0*kl-4 .D0*mul) / (5 .DO* (4 .D0*inul+3 .D0*kl) ) 
sl212= (3 .DO* (2 .D0*mul+kl) ) / (5 .DO* (4 .D0*mul+3 .D0*kl) ) 
sl313= (3 .DO* (2 .D0*mul+kl) ) / (5 .DO* (4 .D0*mul+3 .D0*kl) ) 
s2323=(3.D0* (2.D0*mul+kl) ) /(5.D0* (4 .D0*mul+3 .D0*kl) ) 
go  to  99 
endif 

if  ( ( (a.gt.c) .and. (c.gt.b) ) .or. ( (a.eq.c) .and. (c.gt.b) ) )  then 

c  interchange  2  and  3 

al=a 
bl=b 
cl=c 
b=c 
c=bl 
ip=l 
endif 

if  ( ( (b.gt.a) .and. (a.gt.c) ) .or .( (a.eq.c) .and. (c.lt.b) ) )  then 
c  interchange  1  and  2 

al=a 
bl=b 
cl=c 
a=b 
b=al 
ip=2 
endif 

if  (( (b.gt.c) .and. (c.gt. a) ) .or. ( (b.eq.c) .and. (b.gt.a) ) )  tlien 
c  make  2->l,  3->2,  l->3 

al=a 
bl=b 
cl=c 
a=b 
b=c 
c=al 
ip=3 
endif 

if  ( (c.gt .b) .and. (b.gt.a) )  then 
c  interchange  1  and  3 

al=a 
bl=b 
cl=c 
a=c 
c=al 
ip=4 
endif 

if  ( ( (c.gt.a) .and. (a.gt.b) ) .or. ( (a.eq.b) .and. (b.lt.c) ) )  then 
make  3->l/  l->2,  2->3 
al=a 
bl=b 
cl=c 
a=c 


c 


b=al 

c=bl 

ip=5 

endif 

if  ( (a.eq.b) .and. (b.gt.c) )  then 
call  flush (15) 

ia=2 .DO*p*a*a*c/ (a*a-c*c) 
ia=ia/dsgrt {a*a-c*c) 

ia=ia* (dacos (c/a) -( (c/a) *dsqrt (l.DO- ( (c/a) **2) ) )  ) 
ib=ia 

ic=4.D0*p-ia-ib 

ibc=(ic-ib) / (3 .DO* (b*b~c*c) ) 

iac=(ic-ia) / (3 .DO* (a*a-c*c) ) 

iab=0.25*( ( (4.D0*p) / (3.D0*a*a) ) -iac) 

iba=iab 

ica=iac 

icb=ibc 

iaa=3 .D0*iab 

ibb=( (4.D0*p) / (3.D0*b*b) ) -iba^ibc 
icc=( (4.D0*p) / (3.D0*c*c) )-ica-icb 
go  to  98 
endi  f 

if  ( (b.eq.c) .and. (a.gt.b) )  then 
call  flush (15) 

ib=2 .D0*p*a*c*c/ (a*a-c*c) 
ib=ib/dsqrt (a*a-c*c) 

ib=ib* ( ( (a/c) *dsqrt ( ( (a/c) **2) -l.DO) ) - 
&  dlog( (a/c) +dsqrt ( (a/c) **2-1 .DO)  )  ) 

ic=ib 

ia=4 .D0*p-ib-ic 

iab=(ib-ia) / (3. DO* (a*a-b*b) ) 

iac= (ic-ia) / (3 .DO* (a*a-c*c) ) 

ibc=0 . 25* ( ( (4 .D0*p) / (3 .D0*b*b) ) -iab) 

iba=iab 

ica=iac 

ich-ihc 

iaa= ( (4 .D0*p) / (3 .D0*a*a) ) -iab-iac 
icc=( (4.D0*p) / (3 .D0*c*c) ) -ica-icb 
ibb=3.D0*ibc 
go  to  98 
endif 

c - input  value  of  a;b  and  c 

if  ( (a.gt.b) .and. (b.gt.c) )  then 
ia=4.D0*p*a*b*c/ (a*a-b*b) 
ia=ia/dsqrt ( (a*a-c*c) ) 
x=dsqrt ( ( (a/c) **2) -1 .DO) 
qqc=b*b-c*c 
qqc=qqc/ (a*a-c*c) 
qqc=dsqrt (qqc) 
aa=l .DO 
bb=l.D0 

f f=el2 (x,qqc,aa,bb) 

qqc2=qqc**2 

aa=l .DO 

bb=l.D0 

e=el2 (x,qqc,aa,qqc2) 
ia=ia* (f f-e) 

ic=4 .D0*p*a*b*c/ (b*b-c*c) 

ic=ic/ (dsqrt (a*a-c*c) ) 

ic=ic* ( (b* (dsqrt (a*a-c*c) ) / (a*c) ) -e) 

ib=4 .D0*p-ia-ic 

iab= (ib-ia) / (3 .DO* (a*a-b*b) ) 

ibc= (ic-ib) / (3 .DO* (b*b-c*c) ) 

iac='(ic-ia)  /  (3  .DO*  (a*a-c*c)  ) 

iba=iab 

ica=iac 

icb=ibc 

iaa=( (4.D0*p) / (3.D0*a*a) ) -iab-iac 
ibb= ( (4 .D0*p) / (3 .D0*b*b) ) -iba-ibc 
icc= ( (4 .D0*p) / (3 .D0*c*c) ) -ica-icb 
endif 

c - elements  of  S 

98  q=3.D0/(8.D0*p) 

q=q/ (l.DO-nu) 
r=l.D0-2.D0*nu 
r=r/ (8. D0*p* (l.DO-nu)) 
sllll=q*a*a*iaa+r*ia 
s2222=q*b*b*ibb+r*ib 
s3333=q*c*c*icc+r*ic 
sll22=q*b*b*iab-r*ia 
sll33=q*c*c*iac-r*ia 
s2233=q*c*c*ibc-r*ib 
.  s2211=q*a*a*iba-r*ib 
s3311=q*a*a*ica-r*ic 


s3322=q*b*b*icb-r*ic 

sl212=0.5*q* (a*a+b*b) *iab+0.5*r* (ia+ib) 
sl313=0.5*q* (a*a+c*c) *iac+0.5*r* (ia+ic) 
s2323=0 . 5*q* (b*b+c*c) *ibc+0 . 5*r* (ib+ic) 

99  continue 

if  { spath . eq . . true . )  then 
s{l,l)=sllll 
s(l,2)=sll22 
s(l,3)=sli33 
s(2,l)=:s2211 
s(2,2)=s2222 
s(2,3)=s2233 
s(3,l)=s3311 
s{3,2)=s3322 
s(3,3)=s3333 
s(4,4)=2.D0*sl212 
s(5,5)=2.D0*s2323 
s(6,6)=2.D0*sl313 
do  100  i=l,6 
do  110  j=l,6 
sl(i, j)=s(i, j) 

110  continue 

100  continue 

if  (ip.eq.l)  then 
s(l,2)=sl(l,3) 
s(l,3)=sl(l,2) 
s(2,l)=sl(3,l) 
s(2,2)=sl(3,3) 
s(2,3)=sl(3,2) 
s{3,l)=sl(2,l) 
s(3,2)=sl(2,3) 
s(3,3)=sl(2,2) 
s(4,4)=sl{6,6) 
s(6,6)=sl(4,4) 
else  if  (ip.eq.2)  then 
s(l,l)=sl(2,2) 
s(l,2)=sl(2,l) 
s(l,3)=sl(2,3) 
s(2,l)=sl{l,2) 
s(2,2)=sl(l,l) 
s(2,3)=sl(l,3) 
s(3,l)=sl(3,2) 
s(3,2)=sl(3,l) 
s(5,5)=sl(6,6) 
s(6,6)=sl(5,5) 
else  if  (ip.eq.3)  then 
s(l,l)=sl(3,3) 
s(l,2)=sl(3,l) 
s(l,3)=sl{3,2) 
s(2,l)=sl(l,3) 
s(2,2)=sl{l,l) 
s(2,3)=sl{l,2) 
s(3,l)=sl(2,3) 
s(3,2)=sl(2,l) 
s{3,3)=sl(2,2) 
s(4,4)=sl(6,6) 
s{5,5)=sl(4,4) 
s(6,6)=sl(5,5) 
else  if  {ip.eq.4)  then 
s(l,l)=sl(3,3) 
s(l,2)=sl{3,2) 
s(l,3)=:sl(3,l) 
s(2,l)=sl{2,3) 
s(2,3)=sl(2,l) 
s(3,l)=sl(l,3) 
s(3,2)=sl(l,2) 
s(3,3)=sl(l,l) 
s(4,4)=sl(5,5) 
s(5,5)=sl(4,4) 
else  if  (ip.eq.5)  then 
s(l,l)=sl(2,2) 
s(l,2)=sl{2,3) 
s(l,3)=sl(2,l) 
s{2,l)=sl(3,2) 
s(2,2)=sl(3,3) 
s(2,3)=sl(3,l) 
s(3,l)=sl(l,2) 
s(3,2)=sl(l,3) 
s(3,3)=sl(l,l) 
s(4,4)=sl(5,5) 
s(5,5)=sl(6,6) 
s(6,6)=sl(4,4) 
endif 
else 


c 


elements  of  Pi 
q=3.D0/(8.D0*p) 
q=q/ (l.DO-nu) 
r=l.D0-2 .D0*nu 
r=r/ (8.D0*p* (l.DO-nu) ) 

P1313=(ic-ia) / (8.D0*p) 

P1331=(ic-ia) / (8.D0*p) 

P3131=(ia-ic) / (8.D0*p) 

P3113=(ia-ic) / (8.D0*p) 

P1212= (ib-ia) / (8 .D0*p) 

P1221=(ib-ia) / (8.D0*p) 

P2121= (ia-ib) / (8 .D0*p) 

P2112= (ia-ib) / (8 .D0*p) 

P3232= (ib-ic) / (8 .D0*p) 

P3223=(ib-ic)/(8.D0*p) 

P2323=(ic-ib)/(8.D0*p) 

P2332=(ic-ib)/(8.D0*p) 
ppl212=pl212 

ppl313=pl313 
pp2323=p2323 

ppl221=-ppl212 
pp2121=ppl221 
pp2112=-pp2121 

ppl331=-ppl313 
pp3131=ppl331 
pp3113=-pp3131 
pp2332=-pp2323 
pp3232=pp2332 
pp3223=-pp3232 
if  (ip.eq.l)  then 
pl212=ppl313 
pl313=ppl212 
p2323=pp3232 
else  if  (ip.eq.2)  then 
pl212=pp2121 
pl313=pp2323 
p2323=ppl313 
else  if  (ip.eq.3)  then 
pl212=pp3131 
pl313=pp3232 
p2323=ppl212 
else  if  (ip.eq.4)  then 
pl212=pp3232 
pl313=pp3131 
p2323=pp2121 
else  if  (ip.eq.5)  then 
pl212=pp2323 
pl313=pp2121 
p2323=pp3131 
endif 
endif 

if ( (ip.eq.l) .or. (ip.eq.2) -or. (ip.eq.3) -or. (ip.eq.4) .or. (ip.eq.5) ) 
&  then 
a=al 
b=bl 
c=cl 
endif 
return 
end 

double  precision  FUNCTION  EL2  (X,QQC,  AA.,BB) 
c  implicit  real*8  (a-h,o-z) 

real*8  PI,ca,cb,x,qqc,aa,bb,qc,a,b,c,d,p 
real*8  z,eye,y,ff,em,e,g 
integer  1 

PARAMETER(PI=3. 14159265358979,  CA=1.0D-6,  CB=1.0D-12) 
c  PARAMETER (PI=3. 14159265,  CA=1.0D-6,  CB=1.0D-12) 

IF(X.EQ.0.0)THEN 

EL2=0. 

ELSE  IF(QQC.NE.0.0)THEN 
QC=QQC 
A=AA 
B=BB 
C=X**2 
D=1.D0+C 

P=DSQRT( (l.D0+QC**2*C) /D) 

D=X/D 
C=D/ (2.*P) 

Z=A-B 

EYE=A 

A=0.5*(B+A) 

Y=DABS(1.D0/X) 

FF=0  . 

L=0 


1 


EM=1.D0 

QC=DABS(QC) 

B=EYE*QC+B 

E=EM*QC 

G=E/P 

D=FF*G+D 

FF=C 

EYE=A 

P=G+P 

C=0.5*(D/P+C) 

G=EM 

EM=QC+EM 

A=0.5*(B/EM+A) 

Y=-E/Y+Y 

IF(Y.EQ.O. )Y=DSQRT(E) *CB 
IF(DABS(G-QC) .GT .CA*G) THEN 
QC=DSQRT(E)^2. 

L=L+L 

IF(Y.LT.0.)L=:L+1 
GO  TO  1 
ENDIF 

IF(Y.LT.0.)L=L+1 
E=(DATAN{EM/Y)+PI*L) *A/EM 
IF{X.LT.0.)E=-E 
EL2=E+C*Z 
ELSE 

PAUSE  'failure  in  EL2 ' 
ENDIF 
RETURN 
END 


C  This  is  a  routine  to  invert  matrices  (6x6) 
subroutine  inverse (a, y) 
realms  a(6,6) ,y(6,6) ,d 
integer  indx(6),i,j 
do  10  i=l,6 
do  11  j=1.6 
if  (i.eq.j)  then 
y{i, j)=1.0 
else 

y(i,j)=0.0 

endif 

11  continue 

10  continue 

call  ludcmp (a, 6, 6, indx, d) 
do  13  j=l,6 

call  lubksb{a,6,6,indx,y(l, j) ) 

13  continue 

return 
end 


C  This  is  a  routine  to  invert  matrices  (3x3) 
subroutine  inverse3x3 (a,y) 
real*8  a(3,3) ,y{3,3) ,d 
integer  indx(3),i,j 
do  10  i=l,3 
do  11  j=l,3 
if  (i.eq.j)  then 
y{i, j)=1.0 
else 

y(i,j)=0.0 

endif 

11  continue 

10  continue 

call  ludcmp(a, 3, 3, indx,d) 
do  13  j=l,3 

call  lubksb(a, 3,3,indx,y(l, j ) ) 

13  continue 

return 
end 


SUBROUTINE  LUBKSB  ( A,  N,  NP,  INDX,  B) 
implicit  real*8  (a*-h,  o-^z) 
implicit  integer  (i-n) 
real*8  A(NP,NP) ,B(N) 
integer  indx(n) 

11=0 

DO  12  1=1, N 


LL=INDX(I) 

SUM=B{LL) 

B(LL)=B{I) 

IF  (II.NE.O)THEN 
DO  11  J=II,I-1 

SUM=SUM-A(I,  J)  *B(J) 

11  CONTINUE 

ELSE  IF  (SUM.NE.O.)  THEN 
II=I 
ENDIF 
B(I)=SUM 

12  CONTINUE 

DO  14  I=N,1,-1 
SUM=B(I) 

IF{I.LT.N)THEN 
DO  13  J=I+1,N 

SUM=SUM-A(I,  J)  *B(J) 

13  CONTINUE 
ENDIF 

B(I)=SUM/A{I,I) 

14  CONTINUE 
RETURN 
END 


SUBROUTINE  LUDCMP  ( A,  N,  NP,  I^X,  D)  ^ 
integer  indx(n)  , nitiax,np,n,  i,  j  , k,  imax 
real*8  tiny,d 

PARAMETER  (NMAX=100, TINY=1 . OE-20) 
real*8  A(NP,NP)  ,  W(NMAX)  ,  aaitiax,  sum,duin 
D=l. 

DO  12  1=1, N 
AAMAX=0. 

DO  11  J=1,N 

IF  (DABS(A(I, J) ) .GT.AAMAX)  AAMAX=DABS (A (I , J) ) 

11  CONTINUE 

IF  (AAMAX.EQ.O. )  PAUSE  ' Singular  matrix . ' 
W(I)=1./AAMAX 

12  CONTINXre: 

DO  19  J=1,N 

IF  (J.GT.l)  THEN 
DO  14  1=1, J-1 
SUM=A(I, J) 

IF  (I.GT.l)THEN 
DO  13  K=1,I-1 

SUM=SUM-A(I,K) *A(K, J) 

13  CONTINUE 
A(I, J)=SUM 

ENDIF 

14  CONTINUE 
ENDIF 
AAMAX=0. 

DO  16  I=J,N 
SUM=A(I, J) 

IF  (J.GT.DTHEN 
DO  15  K=1,J“1 

SUM=SUM-A(I,K) *A(K, J) 

15  CONTINUE 
A(I, J)=SUM 

ENDIF 

DUM=W(I)  *DABS(SUM) 

IF  (DUM.GE.AAMAX)  THEN 
IMAX=I 
AAMAX=DUM 
ENDIF 

16  CONTINUE 

IF  (J.NE. IMAX) THEN 
DO  17  K=1,N 
DUM=A{IMAX,K) 

A(IMAX,K)=A(J,K) 

A(J,K)=DUM 

17  CONTINUE 
D=-D 

W{IMAX)=W(J) 

ENDIF 

INDX(J)=IMAX 

IF(J.NE.N)THEN 

IF(A{J,J) .EQ.0.)A(J, J)=TINY 
DUM=1. /A(J, J) 

DO  18  I=J+1,N 

A(I, J)=A(I, J) *DUM 

18  CONTINUE 
ENDIF 

19  CONTINUE 


n  o 


IF{A(N,N)  .EQ.O.)A{N,N)=TINY 

RETURN 

END 


This  is  a  routine  to  multiply  6x6  matrices  (remember  4th  order  tensors 
maybe  written  as  6x6  matrices  without  losing  the  tenosrial  properties) . 
subroutine  tenprod (a, b, c) 
real*8  a(6,6),  b(6,6),  c(6,6) 
integer  i,j,k 
do  20  i=l/6 
do  30  j=l,6 
c(i/ j)=0.0 
30  continue 

20  continue 

do  10  i=l,6 
do  11  j=l,6 
do  12  k=l,6 

c(i, j)=a(i,k) *b{k, j)+c(i,  j) 

12  continue 

11  continue 

10  continue 

return 
end 

subroutine  tenmatprod (m,  s, d) 
real*8  m(6, 6) , s (6) ,d(6) 
integer  i,j 
do  5  i=l,6 
d(i)=0.0 
5  continue 

do  10  i=l,6 
do  20  j=l,6 
d(i)=m(i, j) *s(j)+d{i) 

20  continue 

10  continue 

return 
end 


subroutine  rowcolumnprod(p,q, r) 
real*8  p(6) / q(6) , r 
integer  i 
r=0.0 

*  do  5  i=l,6 

r=r+p(i) *q(i) 

5  continue 

return 
end 


c -  This  is  to  calculate  a  part  of  the  'B'  tensors  for  the  composite 

subroutine  Btensor {a,b, c, ad,bd, cd, c2 ,mul, kl ,mu2 , k2 ,  Bmat) 

realms  Bmat (6, 6) ,L1 (6, 6) ,L2 (6, 6) , Si (6, 6) 

real *  8  Ml ( 6 , 6 ) , mul , mu2 , kl , k2 , nul , nu2 

real*8  a,b, c, ad,bd, cd, cl, c2 

real*8  d (6, 6) , dl (6, 6) 

real*8  pil212 , pi2323 ,pil313 

integer  i,j 

logical' spath, debug 

common  /paths/ spa th 

common  /mdebug/ debug 

nul=0.5* (3 .0*kl”2.0*mul) / (3 . 0*kl+mul) 
nu2=0.5* (3.0*k2-2.0*mu2) / (3.0*k2+mu2) 
cl=1.0-c2 

c  Initialize  the  modulus  tensors 

do  10  i=l,6 
do  20  j=l,6 
Ll(i, j)=0.0 
L2  (i, j)=0.0 
20  continue 

10  continue 

c  Obtain  the  (6x6)  modulus  tensors,  LI  and  L2 

Ll(l,l)=kl+(4.0/3,0) *mul 
Ll(2,2)=kl+(4.0/3.0) *mul 
LI (3,3)=kl+(4.0/3.0) *mul 
Ll(l,2)=kl-(2.0/3.0)*mul 
Ll(2,l)=kl-(2.0/3.0) *mul 
Ll(2,3)=kl-(2.0/3.0) *mul 
Ll{3,2)=kl-(2.0/3.0)*mul 


Ll(3,l)=kl-(2.0/3.0) *mul 
LI (1, 3) =kl- (2 . 0/3 .0) *mul 
Ll(4,4)=2.0*niul 
Ll(5,5)=2.0*mul 
Ll{6,6)=2.0*mul 
L2 (l,l)=k2+{4.0/3.0) *mu2 
L2 (2,2)=k2+(4.0/3,0) *mu2 
L2 (3,3)=k2+ (4. 0/3.0) *mu2 
L2(l,2)=k2-(2.0/3.0) *mu2 
L2 (2 , 1) =k2- (2 . 0/3 . 0) *mu2 
L2 (2 , 3 ) =k2- (2 . 0/3 . 0) *mu2 
L2(3,2)=k2-(2.0/3.0)*inu2 
L2(3,l)=k2-(2.0/3.0)  *inu2 
L2 (l,3)=k2-(2.0/3.0)*mu2 
L2(4,4)=2.0*mu2 
L2(5,5)=2.0*mu2 
L2(6,6)=2.0*mu2 

c  Begin  calculating  B  tensor 

call  inverse (LI, Ml) 
call  tenprod(Ml,L2,d) 
do  30  i=l,6 
d(i,i)=d(i,i)-'1.0 

30  continue 

call  inverse (d,dl) 
do  31  i=l,6 
do  32  j=l,6 
d{i, j)=0.0 
32  continue 

31  continue 

c  calculate  S  tensors 

spa th=. true.  ^  _ _ 

call  s tensor (a, b, c, nul , kl ,inul , Si , pil212 , pi 1313 , pi2323 ) 
do  60  i=l , 6 
do  70  j=l,6 

d(i, j)=Si(i, j)-c2*Si{i, j) 

70  continue 

60  continue 

do  80  i=l,6 
do  90  j=l,6 
d(i, j)-d{i, j)+dl(i,  j) 

90  continue 

80  continue 

call  inverse (d, Bmat) 

return 

end 

doubleprecision  FUNCTION  RTSAFE (FUNCD, XI, X2 ,XACC) 
integer  inaxit,j 

real * 8  xl , x2 , xacc , f 1 , df , f h , xl , xh , swap , dxold , dx , f , temp 

logical  neg, yield, check, cal 

common  /mkdata2 /yield, check, neg 

common  /call /cal 

external  funcd 

PARAMETER  (MAXIT=100) 

CALL  FUNCD (X1,FL,DF) 
if  (neg. eq. .true. )  then 
write (15, *) 'RTSAFEl' 
call  flushdS) 
return 
endif 

CALL  FUNCD (X2,FH,DF) 
if  (neg.eq. . true. )  then 
write (15,*) 'RTSAFE2' 
call  flush (15) 
return 
endif 

IF{FL*FH.GE.O. )  THEN 
c  write (15, *) 'data' ,  xl,x2,fl,fh 

c  call  flush(15) 

PAUSE  'root  must  be  bracketed' 

ENDIF 

IF(FL.LT.0.)THEN 

XL=X1 

XH=X2 

ELSE 

XH=X1 

XL=X2 

SWAP=FL 

FL=FH 

FH=SWAP 

ENDIF 

RTSAFE=.5* (X1+X2) 

DXOLD=DABS ( X2 -Xl ) 


DX=DXOLD 

CALL  FUNCD(RTSAFE,F,DF) 
if  (neg.eq. . true. )  then 
write (15,*) 'RTSAFE2' 
call  flush (15) 
return 
endif 

DO  11  J=1,MAXIT 

IF(( (RTSAFE-XH) *DF-F) * ( (RTSAFE-XL) *DF-F) .GE.O. 
*  .OR.  DABS(2.D0*F) .GT.ABS(DXOLD*DF)  )  THEN 

DXOLD=DX 
DX=0.5*(XH-XL) 

RTSAFE=XL+DX 
IF ( XL . EQ . RTS AFE )  then 
if  (cal. eq. .true. )  then 
write (15,*) 'ITER',J 
endif 
RETURN 
endif 
ELSE 

DXOLD=DX 
DX=F/DF 
TEMP=RTSAFE 
RTSAFE=RTSAFE-DX 
IF ( TEMP . EQ . RTSAFE )  then 
if  (cal .eq. . true. )  then 
write{15,*)  'ITER',J 
endif 
RETURN 
endif 
ENDIF 

IF(DABS(DX) .LT.XACC)  then 
if  (cal .eq. . true. )  then 
write(15,*) 'ITER' , J 
call  flush (15) 
endif 
RETURN 
endif 

CALL  FUNCD (RTSAFE, F,DF) 
if  (neg.eq. .true. )  then 
return 
endif 

IF(F.LT.O.)  THEN 
XL=RTSAFE 
FL=F 
ELSE 

XH=RTSAFE 

FH=F 

ENDIF 

11  CONTINUE 

c  PAUSE  'RTSAFE  exceeding  maximum  iterations' 

write { 15, *) 'RTSAFE  NOT  GOOD' 
call  flush{15) 
neg=.true. 

RETURN 

END 


doubleprecision  function  func(dlam) 
real*8  s(3,3) ,phi,MHS (6, 6) ,Amat(6,6) ,Ce(6,6) 
real*8  sn(3,3),n(3,3) , omega (3, 3) ,sige(3,3) 
real* 8  strainc (6) , sigec (6) , vec (3) ,R(3 , 3 ) 
real* 8  f ,wil,wi2, sci (6) ,dumd(6, 6) ,RT(3, 3) 
real *8  fn,  wiln,  wi2n,dlam,  sigyln,dinc(6) 
real*8  sc(6) ,nc(6) ,omegac(6) ,sig_n(6) 
real*8  mul,kl,mu2,k2>  sigyl, sigma_n(6) 

real* 8  dummy (6) , a,b, c, ad, bd, cd, dummyl (3,3) , duinmy2 (3,3) 

real*8  dummyS (6) , duml , dum2 , snl (3,3) , sigyO 

real*8  dep_plas,  eps^lastic,  deps  (3 , 3)  ,  kll ,  eps^plasticn 

real *8  an,bn, cn, adn,bdn, cdn, el212 , e2323 , el313 

real*8  pil212,pi2323,pil313,bmat (6, 6) ,nul 

real*8  mule, kle,mu2e, ek2e,y,MHSn (6 , 6) ,dum5,duml0 

real *8  cp, sp, ct, st, cs, ss, rot (3 , 3) ,rott(3, 3) ,hrd 

real *8  athetan, apsin, aphin, atheta, aphi, apsi 

real*8  cpn, spn, ctn, stn, csn, ssn, rotn (3 , 3 ) , rotnt (3,3) 

real * 8  loadp ( 6 ) , loading 

real*8  zero (6) , totstrain(6) 

real* 8  rotj (3,3) ,rotjn(3,3) , rota (3, 3) 

integer  i,j,ninc 


logical  path, yield, check, neg, spath,yc, debug 
logical  evolf , evolwi, evolti 
common  /controls/evolf, evolwi, evolti 
common  /nOonodulus  /mul ,  kl ,  mu2 ,  k2 ,  sigyO ,  kll 

common  /mkdatal/f ,wil, wi2, sc, strainc, eps_plastic, sigyl,hrd 
common  /mkdata2 /yield, check, neg 
common  /mkelas/mule, kle,mu2e, ek2e 

common  /mkorient/cp, sp,ct, st, cs,ss,rot,rott,atheta,aphi,apsi 
common  /mkrot/R,RT, rot j 
common  /paths/ spa th 
common  /mdebug/debug 
c 

path= . true . 
yield=. true. 
neg= . false , 

nul=0.5* (3.0*kl-2.0*mul) / (3.0*kl+mul) 
if (debug . eq . . true . )  then 
write  { 15 ,  * )  '  fund ' 
call  flush(15) 
endif 
c 

c=1.000 

a=c/wil 

b=c/wi2 

ad=a 

bd=b 

cd=c 

c 

c -  calculate  Meffective 

call  Meffective (a , b, c , ad, bd, cd, f ,  mul , kll , mu2 , k2 , MHS , path) 
do  100  i=l,6 


do  110  j=l,6 

MHS ( i , j ) =3 . 0  *mul*MHS ( i , j ) 

110  continue 

100  continue 
c 

- - 

c -  Estimate  stress 

c -  STEP  1:  Estimate  the  elastic  predictor 


path= . false. 

call  Meffective(a,b,c,ad,bd,cd, f ,mule,kle,mu2e,ek2e,Ce,path) 

path=.true. 

call  tenmatprod(Ce, strainc, sigec) 
do  201  i=l,6 
loadp (i) =sigec (i) 

201  continue 

siged,  l)=sc(l)+sigec(l) 

sige (2,2) =sc(2) +sigec(2) 

sige (3,3) =sc ( 3 ) +sigec { 3 ) 

sige ( 1 , 2 ) = ( sc ( 4 ) +sigec ( 4 ) ) /dsqr t ( 2 . DO ) 

sige (2,3) = (sc (5) +sigec (5) ) /dsqrt (2 .DO) 

sige ( 3 , 1 ) = ( sc ( 6 ) +sigec ( 6 ) ) /dsqrt ( 2 . DO ) 

sige (2, 1) =sige(l,2) 

sige(3,2)=sige(2,3) 

siged, 3)=sige(3,l) 

sigecd)  =sige(l,  1) 

sigec (2) =sige(2 , 2) 

sigec(3) =sige(3, 3) 

sigec (4) =sige (1 , 2) *dsqrt (2 .DO) 

sigec (5) =sige (2 , 3) *dsqrt(2.D0) 

sigec(6)=sige(3, 1) *dsqrt(2.D0) 

c -  STEP  la:  Determining  the  stress  on  the  yield  surface 

sigma_n (1 ) =sc (1) 
sigmai.n  ( 2 )  =sc  ( 2 ) 
sigma_n ( 3 ) =sc ( 3 ) 
sigma_n ( 4 ) =sc ( 4 ) 
sigma_n ( 5 ) =sc ( 5 ) 
sigma_n(6) =sc (6) 

call  yieldtest (MHS, sigma^n, sigyl , f ,yc,y) 
if  (yc.eq. . false . )  then 
call  yieldtest (MHS, sigec, sigyl, f,yc,y) 
if  (yc.eq. .true. )  then 

c  this  is  the  step  where  the  behavior  becomes  plastic  (the  previous 

c  step  was  elastic) 

c  write (15, *) 'Becoming  plastic' 

do  127  i=l,6 
zero(i)=0.0 

totstrain ( i ) =strainc ( i ) 

127  continue 

ninc=100.0 

call  stre (sigiaa_n, sigec, zero, totstrain, sig_n, Ce, nine) 
else 

write (15,*) 'NOT  PLASTIC' 
call  flushdS) 


130 


do  130  i=l,6 
sig_n (i) =sc (i) 
continue 
endif 
else 

call  flushdS) 
do  131  i=l,6 
sig_n(i)=sc(i) 


131  continue 

endif 

c -  calculate  N(3,3) 


call  tenmatprod(MHS, sig_n,nc) 
do  200  i=l,6 
nc(i)=nc(i) / (1.0-f ) 
nc(i)=2 ,0*nc(i) /sigyl 
dummy ( i ) =nc ( i ) 
dummyS  ( i )  =nc'(  i ) 

200  continue 

n(l,l)=nc(l) 

n(2,2)=nc(2) 

n(3,3)=nc(3) 

n(l,2)=nc(4) /dsqrt(2 .DO) 
n(2,l)=n(l,2) 
n(2,3)=nc(5) /dsqrt(2 .DO) 
n(3,2)=n(2,3) 
n(3,l)=nc{6) /dsqrt(2.D0) 
n(l,3)=n(3,l) 

call  rowcoluinnprod(loadp,dummy5,  loading) 
if  { loading. le. 0.0)  then 
write (15, *)' loading  in  func loading 


endif 

do  202  i=:l,6 
dummyS ( i ) =nc ( i ) 

202  continue 

c -  STEP  2:  Remaining  part  of  the  stress 

c - STEP  2a:  Calculating  'omega' 

spath= . false . 


do  17  i=l,6 
do  18  j=l,6 
dumd (i , j ) =0 . 0 
18  continue 

17  continue 

if (debug . eq. . true . )  then 

write (15, *) 'func3' ,a,b,c,nul,kl,mul,dumd,pil212,pil313, 
&  pi2323 

call  flushdS) 
endif 

call  stensor (a, b, c, nul , kl ,mul, dumd, pil212 ,pil313,pi2323) 

el212=pil212-f*pil212 

e2323=pi2323-f*pi2323 

el313=pil313-f*pil313 

spath= . true. 

i f ( debug . eq . , true . )  then 

write (15, *) 'func3_a' , ad, bd, cd, a,b, c 

call  flushdS) 

endif 

call  Btensor (a, b, c, ad, bd, cd, f ,mul , kl , mu2 , k2 , Bmat) 
do  41  i=l,3 
do  42  j=l,6 
Bmatd,  j)=0.0 
42  continue 

41  continue 

Bmat  ( 4 , 1 )  =*-2 . 0 *el212 *Bmat  ( 4 , 1 ) 
Bmat(4-,2)=-2.0*el212*Bmat(4,2) 

Bmat(4,3)=-2.0*el212*Bmat (4,3) 

Bmat (4, 4) =-2 . 0*el212*Bmat (4,4) 

Bmat(4, 5) =-2 .0*el212*Bmat (4, 5) 

Bmat  ( 4 , 6 )  =-2 . 0*el212*Binat  ( 4 , 6 ) 
Bmat(5,l)=-2.0*e2323*Bmat(5,l) 
Bmat(5,2)=-2.0*e2323*Bmat(5,2) 

Bmat ( 5 , 3 ) =-2 . 0*e2323*Bmat ( 5 , 3 ) 
Bmat(5,4)=-2.0*e2323*Bmat(5,4) 

Bmat ( 5 , 5 ) =-2 . 0*e2 323 *Bmat ( 5 , 5 ) 

Bmat (5, 6) =-2 .0*e2323*Bmat (5, 6) 

Bmat(6,l)=-2.0*el313*Bmat (6,1) 

Bmat(6,2)=-2.0*el313*Bmat (6,2) 
Bmat(6,3)=-2.0*el313*Broat(6,3) 
Bmat(6,4)=-2.0*el313*Bmat(6,4) 

Bmat (6, 5) =-2 .0»el313*Bmat (6, 5) 
Bmat(6,6)=-2.0*el313*Bmat(6,6) 
if  (debug. eq.  .true. )  tlien 
writedS,  *)  'func3_b'  ,a,b,c,ad,bd,cd 
call  fluslidS) 
endif 


o  n 


call  Atensor (a,b, c,ad,bd, cd, f ,mul,kl,mu2,k2, Amat) 
c - 

if (debug . eq . . true . )  then 
write (15, *) 'func4' 
call  flush (15) 
endif 

do  3  i=l,3 
do  4  j=l/3, 
omegad,  j)  =0.0 
4  continue 

3  continue 

do  1  i=l,3 

omega(l,2)=-Binat  (4,  i)  *nc  (i) +oinega  (1, 2) 
omega (2,3) =-Bmat ( 5 , i ) *nc ( i ) +omega (2,3) 
omega(l,3)=“Binat(6,i)  *nc  (i) +omega  (1, 3) 

1  continue 
do  2  i=4,6 

omegad,  2)  =-Bmat  (4,  i)  *nc  (i) +omega  (1, 2 ) 
omega(2,3)=-Bmat(5,i) *nc (i) +omega (2 , 3 ) 
omega (1 , 3 ) =-Bmat ( 6 , i) *nc ( i) +omega (1 , 3 ) 

2  continue 

omega (2,1) =“Omega (1,2) 
omega (3,2) =-omega (2,3) 
omega ( 3 / 1 ) =-omega (1,3) 
do  7  i=l,3 
do  8  j=l,3 

omega (i , j ) =omega (i , j ) /dsqr t ( 2 . OdO ) 

8  continue 

7  continue 

call  tenmatprod(Amat,nc,dinc) 
if  (dabs(a-b) .gt.0.01)  then 

omega (1,2) =omega (1,2)- (a*a+b*b) *dinc ( 4 ) / (dsqrt ( 2 . OdO ) * ( a*a-b*b) ) 
endif 

if  (dabs(a-c) .gt.0.01)  then 

omegad,  3)=omega(l,  3)  -  (a*a+c*c)  *dinc  (6)  /  (dsqrt  (2  .OdO)  *  (a*a-c*c)  ) 
endif 

if  (dabs (c-b) .gt. 0.01)  then 

omega(2, 3) =omega(2,3) - (b*b+c*c) *dinc (5) / (dsqrt (2 . OdO) * (b*b-c*c) ) 
endif 

omega (2,1) =-omega (1, 2 ) 

omega (3 , 2 ) =-omega (2 , 3 ) 

omega (3,1) =-omega (1,3) 

omegac ( 1 ) =omega (1,1) 

omegac ( 2 ) =omega (2,2) 

omegac ( 3 ) =omega ( 3 , 3 ) 

omegac ( 4 ) =omega (1,2) *dsqr t ( 2 . DO ) 

omegac (5) =omega(2, 3) *dsqrt (2 .DO) 

omegac (6) =omega (1 , 3 ) *dsqrt (2 . DO) 

c - 

s(l,l)=sc(l) 
s(2,2)=sc(2) 
s(3,3)=sc(3) 

s ( 1 , 2 ) =sc ( 4 ) /dsqrt ( 2 , DO ) 
s(2,3)=sc(5) /dsqrt (2. DO) 
s  ( 3 , 1 ) =sc ( 6 ) /dsqrt ( 2 . DO ) 
s(2,l)=s(l,2) 
s(3,2)=s(2,3) 
s(l,3)=s(3,l) 
write (15, *) 'In  func',sc 
call  flush(15) 

call  matprod(s, omega, dummy 1) 
call  matprodiotnega,  s,  dummy2) 
do  300  i=l,3 
do  310  j=l,3 

sn  (i ,  j )  =duminyl  d ,  j )  -dummy2  d ,  j  ) 

CCC  sn(i,j)=0.0 

310  continue 

300  continue 

do  450  i=l,6 
dummy5 ( i ) =nc ( i ) 

450  continue 

if (debug . eq . . true . )  then 
write ( 15 , * ) ' f unc5 ' 
call  flush (15) 
endif 

call  tenmatprod (Ce, dummy 5, dummy) 
sn(l,l)=sn(l, 1) -dummy (1) 
sn(2,2)=sn(2,2) -dummy(2) 
sn(3, 3) =sn(3, 3) -dummy (3) 
sn(l,2)=sn(l,2) -dummy (4) /dsqrt (2. DO) 
sn(2,3)=sn(2, 3) -dummy(5) /dsqrt (2 .DO) 
sn(3, 1) =sn(3, 1) -dummy (6) /dsqrt (2 .DO) 
sn(2,l)=sn(l,2) 
sn(3,2)=sn{2,3) 


sn(l,3)=sn(3,l) 
do  320  i=l,3 


do  330  j=l,3 
duinmyl  ( i  /  j )  =sn  ( i ,  j  ) 

330  continue 

320  continue 

Q - STEP  3 :  put  them  together 


snl(l,l)=sige(l,l)+dlam*sn(l,l) 
snl (2,2) =sige (2,2) +dlam*  sn ( 2 , 2 ) 
snl (3,3) =sige (3,3) +dlam*  sn ( 3 , 3 ) 
snl (1,2) =sige (1,2) +dlam*sn (1,2) 
snl (2, 3) =sige(2, 3) +dlam*sn (2 , 3 ) 
snl(3,l)=sige(3, 1) +dlam*sn (3 , 1) 
snl(2,l)=snl(l,2) 
snl (3, 2) =snl (2,3) 
snl(l,3)=snl(3,l) 

Q _  These  components  of  snl  are  now  NOT  in  the  coordinate  farme 

c -  coincides  with  the  orientation  of  the  inclusions.  This  is  due 

c _  the  algorithm  used  (see  notes  for  more  details) .  In  order  to 

c  obtain  them  in  this  coordinate  frame: 

call  matprod(R, snl, dummy 1) 
call  matprod ( dummy  1,  RT,  snl) 
c  write (15,*) ' snl ' , snl (3 , 3 ) 

ccall  flush(15) 

- - 7“”7  7 

c -  estimate  the  new  state  variables  for  this  iteration 

if (debug . eq . . true . )  then 
writedS,  *)  'func6' 
call  flush(15) 
endif 

call  Atensor (a, b, c, ad, bd, cd, f ,mul , kl,mu2 , k2 , Amat) 

if (debug . eq. . true . )  then 
write ( 15 , * ) ' func6a ' 
call  flush(15) 
endif 

if  (evolf .eq. . true. )  then 
fn=f+(n(l,l)+n(2,2)+n(3,3) ) *(1.0-f) *dlam 
else 
fn=f 
endif 

if  ( (fn.lt. 0.0) .or. (fn.gt. 1.0) )  then 
write ( 15, *) 'porosity  estimate  wrong' , f , dlam 
call  flush(15) 
c  neg=.true. 

c  return 

f=0.0 
endif 

dummy  ( 1 )  =Amat  (3,1)  -Amat  (1,1) 
dummy { 2 ) =Amat ( 3 , 2 ) -Amat ( 1 , 2 ) 
dummy ( 3 ) =Amat (3,3) -Amat (1,3) 
diammy  ( 4 )  =Amat  (3,4)  -Amat  (1,4) 
dummy  ( 5 )  =Amat  (3,5)  -Amat  (1,5) 
dummy  ( 6 )  =Amat  (3,6)  -Amat  (1,6) 
call  rowcolumnprod  (dummy ,  duinmy5 ,  duml ) 
do  400  i=l,6 
dummy5 ( i ) =nc ( i ) 

400  continue 

dummy ( 1 ) =Amat (3,1) -Amat (2,1) 
dummy ( 2 ) =Amat (3,2) -Amat (2,2) 
dummy ( 3 ) =Amat (3,3) -Amat (2,3) 
dummy (4 ) =Amat (3 , 4 ) -Amat (2,4) 
dummy ( 5 ) =Amat (3,5) -Amat (2,5) 
dummy ( 6 ) =Amat (3,6) -Amat (2,6) 
call  rowcolumnprod  ( d\immy ,  dummy 5 ,  dum2 ) 
if (debug. eq. . true. )  then 
write { 15 , * ) ' f unc6b ' 
call  flush (15) 
endif 

'if  (evolwi .eq. . true. )  then 
wi  ln=wi  1 + dl  am*  wi  1  *  duml 
wi2n=wi2+dlam*wi2*dum2 
else 

wiln=wil 

wi2n=wi2 

endif 

dum5=duml 

duml0=dxim2 

if  ( (wiln.le.0.0) .or. (wi2n. le . 0 . 0) )  then 
write (15, *) 'aspect  estimate  wrong' , 'dlam=' , dlam 
call  flush (15) 
neg= . true . 
return 
endif 

estimating  new  orientations 


c 


if  (evolti.eq. . true. )  then 
athetan=atheta- (omega (2, 3) *cp+omega (1, 3) *sp) *dlam 
if  (dabs (wi2n-1.0) .It. 0.01)  then 
athetan=0 .0 
endif 

if  (st.ne.0.0)  then 

aphin=aphi- ( (omega (2, 3) *ss/st) -omega (1, 3) *cs/st) *dlam 
apsin=apsi+ (-omega (1, 2) + (ct/st )* (omega (2,3) *ss-omega( 1,3) *cs) ) * 

Sc  dlam 

else 

apsin=0 . 0 
aphin=0 . 0 
endif 

call  matprod(R, rotj , rotjn) 
else 

aphin=:aphi 
apsin=apsi 
athetan=atheta 
do  27  i=l,3 
do  28  j=l,3 
rotjn(i, j)=rotj (i, j) 

28  continue 

27  continue 

endif 

ctn=dcos (athetan) 

stn=dsin (athetan) 

cpn=dcos (aphin) 

csn=dcos (apsin) 

spn=dsin (aphin) 

ssn=dsin (apsin) 

rota (1,1) =csn*cpn-ssn*ctn*spn 

rota (1,2) =-csn*spn-ssn*ctn*cpn 

rota ( 1 , 3 ) =ssn*stn 

rota (2,1) =ssn*cpn+csn*ctn*spn 

rota (2, 2) =-ssn*spn+csn*ctn*cpn 

rota (2 , 3 ) =-csn*stn 

rota (3,1) =stn*spn 

rota (3,2) =stn*cpn 

rota(3, 3) =ctn 

call  matprod (rota, rot jn, rotn) 
do  77  i=l,3 
do  88  j=l,3 
rotnt ( i , j ) =rotn ( j , i ) 

88  continue 

77  continue 

c - estimating  new  sig_yield 

do  501  i=l,3 
do  502  j=l,3 
deps ( i , j ) =dlam*n ( i ,  j ) 

502  continue 
501  continue 

duml= ( deps (1,1) +deps (2,2) +deps ( 3 , 3 ) ) / 3 . 0 
do  503  i=l,3 

deps ( i , i ) =deps ( i , i ) -duml 

503  continue 

dep_plas=deps (1,1)*  *2  +deps ( 2 , 2 ) *  * 2  +deps ( 3 , 3 ) *  *2 
dep jlas=dep_plas+2 . 0  * (deps (1,2) *  *2+deps (2,3) *  *2  +deps ( 3 , 1 ) *  *2 ) 
dep_plas= (2 . 0/3 . 0 ) *dep_plas 
eps_plasticn=eps_plastic+dsqrt (dep_plas) 

CCC  sigyln=sigy0* ( (1 . 0+eps_plasticn) ** (1.0/3. 0) ) 

sigyln=sigyl 

c - 

C -  calculate  the  yield  function  with  these  new  values  of  the  stress 

c  and  state  variables..... 

c 

c 

cn=1.000 

an=cn/wiln 

bn=cn/wi2n 

adn=an 

bdn=bn 

cdn=cn 

call  Meffective(an,bn, cn, adn,bdn, cdn, fn,mul,kll,mu2,k2,MHSn,path) 
do  500  i=l,6 
do  510  j=l,6 

MHSn(i,  j)=3.0*mul*MHSn(i,  j) 

510  continue 
500  continue 

c -  The  stress  must  now  be  expressed  w.r.t  the  new  basis  —  the  basis 

c  which  is  aligned  with  the  current  (new)  orientation  of  the  voids.  This  is 

c  the  basis  in  which  the  components  of  MHSn  are  obtained, 

c 

c - ;  This  expresses  snl  in  GLOBAL  frame. 

call  matprod (rot, snl, dummy) 


c 


call  matprod (dummy, rot t,snl) 

This  expresses  snl  w.r.t  new  orientation  of  particles 
call  matprod(rotnt, snl, dummy) 
call  matprod (dummy, rotn, snl) 

sci (1) =snl (1,1) 

sci (2) =snl (2,2) 

sci P) =snl (3, 3) 

sci (4) =snl (1, 2) *dsqrt (2 .DO) 

sci (5) =snl (2, 3) *dsqrt (2 .DO) 

sci (6) =snl (3,1) *dsqrt (2 .DO) 

call  tenmatprod(MHSn, sci, dummy) ^ 

call  rowcoluinnprod( sci,  dummy, phi) 

phi=phi/( (1.0-fn) *sigyln) 

phi=phi”Sigyln 

func=phi 

return 

end 


doubleprecision  FUNCTION  RTNEWT (FUNCD, X2 , XACC) 
integer  jmax, j 

c  implicit  doubleprecision (a-h, o^z) 

real*8  x2,xacc, f ,df ,dx 
logical  neg, check, yield 
common  /mkdata2 /yield, check, neg 
external  funcd 
PARAMETER  (JMAX=20) 
rtnewt=x2 
DO  11  J=1,JMAX 

CALL  FUNCD (RTNEWT, F,DF) 
if  (neg. eq. .true. )  then 
return 
endif 
DX=F/DF 

RTNEWT=RTNEWT-DX 

c  RTNEWT=RTNEWT- 0 .000001 

IF ( DABS ( DX ) . LT . XACC )  then 
RETURN 
endif 

11  CONTINUE 

c  PAUSE  'RTNEWT  exceeding  maximum  iterations' 

write (15, *) 'RTNEWT  NOT  GOOD' 
call  flush (15) 
neg=.true. 

END 


subroutine  guessmaker (dlam) 
real*8  MHS (6, 6) ,Ce (6, 6) 
real*8  n(3, 3) ,sige(3, 3) ,L 
real*8  strainc (6) , sigec (6) 
real*8  f,wil,wi2 

real*8  cp, sp, ct, st, cs, ss, rot (3 , 3) ,rott(3,3) , atheta, aphi, apsi 

real*8  dlam, sigyO 

real* 8  sc (6) , nc (6) , sig_n(6) 

real * 8  mul , kl , mu2 , k2 , sigyl , sigma_n ( 6 ) 

real*8  dummy(6) ,a,b,c,ad,bd,cd 

real*8  dummyS (6) , kll,mule, kle,mu2e, ek2e, eps^plastic 
real* 8  phil , f tes t , phi2 , dphidf , h, wiltes t , ate , bte 
real * 8  dphidwi 1 , wi2 test , dphidwi2 , dum5 , duml 0 , y 
real *  8  omega (3,3) , omegac ( 6 ) , dummyl (3,3) , dummy2 (3,3) 
real *8'  bmat (6, 6) ,amat (6, 6) ,dumd(6, 6) , dine (6) ,hrd 
real*8  pil212,pi2323,pil313,el212,e2323,el313 
real* 8  nul, zero(6) , totstrain(6) 
integer  i,j,ninc 

logical  path, yield, check, neg, yc, spath, debug 
logical  evolf , evolwi, evolti 
common  /controls/evolf, evolwi , evolti 
common  /mkmodulus/mul , kl , mu2 , k2 , sigyO , kll 
common  /mkelas/mule, kle,mu2e, ek2e 

common  /mkdatal/f, wil ,wi2 , sc, strainc, eps_plastic, sigyl ,hrd 
common  /mkorient/cp, sp, ct, st, cs, ss, rot , rott, atheta, aphi, apsi 
common  /mkdata2 /yield, check, neg 
common  /mdebug/ debug 


c 

path= . true . 

yield=. true. 

c=1.000 

a=c/wil 

b=c/wi2 

ad=a 


bd=b 

cd=c 


c 

c - 

c -  Estimate  stress 

c - STEP  1:  Estimate  the  elastic  predictor 


path=. false. 

call  Mef f ective (a, b, c , ad, bd, cd, f , mule, kle, mu2e, ek2e, Ce , path) 
path= . true . 

call  tenmatprod(Ce, strainc, sigec) 
siged,  1)  =sc(l)+sigec(l) 
sige (2,2) =sc  d) +sigec{2) 
sige  d  r  3 ) =sc { 3 ) +sigec ( 3 ) 
sige (1, 2) = (sc (4) +sigec (4) ) /dsqrt (2 .DO) 
sige(2, 3)=(sc(5) +sigec(5) ) /dsqrt (2. DO) 
sige ( 3 , 1 ) = ( sc ( 6 ) +sigec ( 6 ) ) /dsqrt (2 . DO ) 
sige (2,1) =sige (1, 2) 
sige ( 3 , 2 ) =s ige ( 2 , 3 ) 
sige (1,3) =sige (3,1) 
sigec(l)=sige(l,l) 
sigec(2) =sige (2, 2) 
sigec(3) =sige (3, 3) 
sigec (4) =sige (1, 2) *dsqrt (2 .DO) 
sigec ( 5 ) =sige ( 2 , 3 ) *dsqr t ( 2 . DO ) 
sigec (6) =sige (3, 1) *dsqrt (2 .DO) 
c - Test  for  yielding... 

Q -  calculate  yield  function  with  stress=elastic  predictor 

call  Mef f ective (a, b, c, ad, bd, cd, f ,mul , kl ,mu2 , k2 ,MHS, path) 
do  101  i=l,6 
do  111  j=l,6 

MRS ( i , j ) =3 . 0  *mul *MHS ( i , j ) 

111  continue 

101  continue 

call  yieldtest(MHS, sigec, sigyl, f,yc,y) 
if  (yc.eq. . false . )  then 
yield= . false . 
go  to  999 
endif 

if (debug. eq. . true. )  then 
write (15, * ) 'guess2 ' 
call  flushdS) 
endif 

c -  STEP  la:  Determining  the  stress  on  the  yield  surface 

sigina__n  (1 )  =sc  (1 ) 
sigma_n (2) =sc (2) 
sigma_n(3) =sc (3) 
sigma_n (4) =sc (4) 
sigTna_n  ( 5 )  =sc  ( 5 ) 
sigma_n(6) =sc (6) 

call  yieldtest (MHS, sigma_n, sigyl, f,yc,y) 
if  (yc . eq .. false . )  then 
cal 1  yieldtest (MHS , sigec , sigyl , f , yc , y ) 
if  (yc . eq. . true . )  then 

c  this  is  the  step  where  the  behavior  becomes  plastic  (the  previous 

c  step  was  elastic) 

do  127  i=l,6 
zero (i) =0 . 0 

totstrain(i) =strainc (i) 

127  continue 

ninc=100.0 

call  stre {sigma_n, sigec, zero, totstrain, sig__n, Ce, nine) 
else 

do  130  i=l,6 
sig_n(i) =sc (i) 

130  continue 
endif 

else 

do  131  i=l,6 
sig_n (i) =sc (i) 

131  continue 
endi  f 

Q -  calculate  N(3,3) 

call  tenmatprod(MHS, sig_n,nc) 
do  200  i=l,6 
nc (i) =nc (i) / (1 . 0-f ) 
nc ( i ) =2 . 0*nc ( i ) /sigyl 
200  continue 

n(l,l)=nc(l) 

n(2,2)=nc(2) 

n(3,3)=nc(3) 

n(l,2)=nc(4) /dsqrt(2.D0) 
n(2,l)=n(l,2) 
n(2,3)=nc(5) /dsqrt (2. DO) 
n(3,2)=n(2,3) 


n  (3, 1)  =nc  (6)  /dsqrt  (2  .DO) 
n(l,3)=n(3,l) 

- - 

c - calculate  guess 

do  410  i=l,6 
duitimyS  ( i )  =nc  ( i ) 

410  continue 

call  tenmatprod(Ce,duinniy5, dummy) 
do  420  i=l,‘6 
dummy 5 ( i ) =nc ( i ) 

420  continue 

call  rowcolumnprod  (dummy,  dimnmyS ,  L ) 
i f ( debug . eq . . true . )  then 
writedS,  *)  'guess3' 
call  flushdS) 
endif 

c -  calculating  omega 

spath= . false . 

nul=0.5*(3.0*kl-2.0*mul) / (3.0*kl+mul) 

call  stensor (a,b,c,nul,kl,mul,dumd,pil212,pil313,pi2323) 

el212=pil212-f*pil212 

e2323=pi2323-f*pi2323 

el313=pil313-f*pil313 

spath=. true. 

call  Btensor (a, b, c, ad, bd, cd, f , mul , kl , mu2 , k2 , Bmat) 
do  41  i=l,3 
do  42  j=l,6 
Bmat (i, j ) =0 . 0 
42  continue 

41  continue 

Bmat (4, 1) =2 . 0*el212*Bmat (4, 1) 

Bmat (4,2) =2.0*el212*Bmat (4,2) 

Bmat  ( 4 , 3 )  =2 . 0 *el2 12 *Binat  ( 4 , 3 ) 

Bmat (4, 4)  =2 .0*el212*Bmat (4, 4) 

Bmat (4,5)=2.0*el212  *Bmat (4,5) 

Bmat ( 4 , 6 ) =2 . 0 *el2 12 *Bmat ( 4 , 6 ) 

Bma t ( 5 , 1 ) =2 . 0 * e2 3 2 3 *Bmat ( 5 , 1 ) 

Bmat  ( 5 , 2 )  =2 . 0*e2323 *Binat  ( 5 , 2 ) 

Bmat (5, 3) =2 . 0*e2323*Bmat (5, 3) 

Bmat(5,4)=2.0*e2323*Bmat(5,4) 

Bmat(5, 5)=2.0*e2323*Bmat (5,5) 

Bmat (5,6)=2.0*e2323*Bmat (5, 6) 

Bmat ( 6 , 1 ) =2 , 0 *el313 *Bmat ( 6 , i ) 

Bmat (6, 2) =2 . 0*el313*Bmat ( 6 , 2 ) 

Bmat ( 6 , 3 ) =2 . 0*el3 13  *Bmat (6,3) 

Bmat (6, 4) =2 .0*el313*Bmat (6, 4) 

Bmat(6, 5) =2.0*el313*Bmat(6,5) 

Bmat (6,6) =2.0*el313*Bmat(6, 6) 

call  Atensor (a,b, c, ad,bd, cd, f ,mul,kl,mu2, k2, Amat) 
do  3  i=l,3 
do  4  j=l,3 
omega (i , j ) =0 . 0 
4  continue 

3  continue 

do  1  i=l,3 

omega (1,2) =-Bmat ( 4 , i ) *nc ( i ) +omega (1,2) 
omega ( 2 , 3 ) =-Bmat ( 5 , i ) *nc ( i ) +omega ( 2 , 3 ) 
omega (1,3)= -Bmat { 6 , i ) *nc ( i ) +omega (1,3) 

1  continue 
do  2  i=4,6 

omega (1,2) =-Bmat ( 4 , i ) *nc ( i ) +omega (1,2) 
omega (2,3) =-Bmat ( 5 , i ) *nc ( i ) +omega (2,3) 
omega(l,3)=-Bmat (6, i) *nc(i)+omega(l,  3) 

2  continue 

omega (2,1) =-omega (1,2) 
omega (3,2) =-omega (2,3) 
omega (3,1) =-omega (1,3) 
do  7  i=l,3 
do  8  j=l,3 

omega ( i , j ) =omega ( i , j ) /dsqrt (2 . OdO ) 

8  continue 

7  continue 

call  tenmatprod (Amat, nc, dine) 
if  (dabs (a-b) .gt. 0.01)  then 

omega(l,2)=omega(l,2) - (a*a+b*b) *dinc(4) / (dsqrt (2 . OdO) * (a*a-b*b) ) 
endif 

if  (dabs(a-c) .gt.0.01)  then 

omega(l,3)=omega(l, 3) -(a*a+c*c) ♦dine (6) / (dsqrt (2 . OdO) * (a*a-c*c) ) 
endif 

if  (dabs (c-b) .gt . 0 . 01)  then 

omega(2, 3)=omega(2, 3) - (b*b+c*c) *dinc(5) / (dsqrt (2 . OdO) * (b*b-c*c) ) 
endif 

omega (2,1) =-omega (1,2) 
omega (3,2) =-omega (2,3) 


omega (3,1) = -omega (1,3) 
omegac ( 1 ) =omega (1,1) 
omegac ( 2 ) =omega (2,2) 
omegac  ( 3 )  =omega  (3,3) 
omegac ( 4 ) =omega (1,2) *dsqr t ( 2 . DO ) 
omegac ( 5 ) =omega (2,3) *dsqr t ( 2 . DO ) 
omegac (6) =omega(l, 3) *dsqrt (2 .DO) 
call  matprod(sig_n,  omega,  dummy  1) 
call  matprod (omega,  sig_n,  durnmy2 ) 
do  12  i=l,3 
do  13  j=l,3 

diommyl  ( i ,  j  )  =diimmyl  ( i ,  j  )  -dummy2  ( i ,  j ) 

13  continue 

12  continue 

do  14,  i=l,3 
do  15  j=l,3 
L=L-dummyl ( i , j ) *n ( i , j ) 

15  continue 

14  continue 

c -  calculating  H 

i f ( debug . eq . . true . )  then 
write (15, *) 'guess4' 
call  flush(15) 
endif 

c - step  A  —  calculating  dphi/df 

H=0.0 

if  (evolf .eq. . true. )  then 
path=. true. 

call  Mef fective(a,b, c,ad,bd,cd, f ,mul, kll,mu2, k2,MHS,path) 
do  1001  i=l,6 
do  1002  j=l,6 
MHS{i, j ) =3 . 0*mul*MHS(i, j ) 

1002  continue 

1001  continue 

call  yieldtest (MHS , sc , sigyl , f , yc , phil ) 
ftest=f+0 . 001*f 

call  Mef fective(a,b, c,ad,bd, cd, f test, mul. kll,mu2,k2, MHS, path) 
do  1003  i=l,6 
do  1004  j=l,6 
MHS ( i , j ) =3 . 0*mul*MHS (i , j ) 

1004  continue 

1003  continue 

call  yieldtest(MHS, sc, sigyl, f test, yc,phi2) 
dphidf= (phi2“phil) / (ftest-f ) 

H=0.0 

H=-dphidf * (1.0-f ) * (nc(l)+nc(2)+nc(3) ) 
endif 

c -  step  B  —  calculating  dphi/dwil 

if  (evolwi .eq. . true. )  then 
wiltest=wil+0 . 001*wil 
ate=c/wiltest 

call  Meffective(ate,b,c,ate,b,c, f ,raul,kll,mu2,k2, 

&  MHS, path) 

do  2003  i=l,6 
do  2004  j=:l,6 
MHS ( i , j ) =3 . 0  *mul ♦MHS ( i , j ) 

2004  continue 

2003  continue 

call  yieldtest(MHS,sc,sigyl, f ,yc,phi2) 
dphidwil= (phi2-phil) / (wiltest-wil) 
dummy ( 1 ) =Amat (3 , 1 ) -Amat (1 , 1) 
dummy ( 2 ) = Ama t ( 3 , 2 ) - Ama t ( 1 , 2 ) 

dummy ( 3 ) =Amat (3,3) -Amat (1,3) 
dummy  ( 4 )  =Ainat  ( 3 , 4 ) -Amat  ( 1 , 4 ) 
dummy  ( 5 )  =Amat  (3,5)  -Amat  (1,5) 
dummy ( 6 ) =Amat ( 3 # 6 ) -Amat ( 1 , 6 ) 
do  400  i=l,6 
duinmy5  ( i )  =nc  ( i ) 

400  continue 

call  r owe olumnprod( dummy, dummy 5, dum5) 

H=H-dphidwil*wil*dum5 

endif 

c -  step  C  —  calculating  dphi/dwi2 

if  (evolwi . eq. . true. )  then 
wi2test=wi2+0 . 001*wi2 
bte=c/wi2test 

call  Mef fective (a,bte, c,a,bte,c, f ,mul,kll,mu2,k2, 

&  MHS, path) 

do  3003  i=l,6 
do  3004  j=l,6 
MHS (i, j ) =3 . 0*mul*MHS (i , j ) 

3004  continue 

3003  continue 

call  yieldtest (MHS, sc, sigyl , f ,yc,phi2) 


o  o 


dphidwi2=  (phi2-phil)  /  (wi2test-wi2) 
do  401  i=l, 6 
dummy 5 ( i ) =nc ( i ) 

401  continue 

dummy  ( 1 )  =Ainat  (3,1)  -"Amat  (2,1) 
dummy  ( 2 )  =Ainat  (3,2)  -  Amat  (2,2) 
dummy  ( 3 )  =Amat  (3,3)  -Amat  (2,3) 
dummy  ( 4 )  =Airiat  (3,4)  -Amat  (2,4) 
dummy  ( 5 )  =Amat  (3,5)  -Amat  (2,5) 
dummy  ( 6 )  =Ainat  (3,6)  -Amat  (2,6) 
call  rovcolumnprod (dummy ,  dummy 5 ,  dumlO ) 
H=H-dphidwi2  *wi2  *duml0 
endif 


L=L+H 


do  450  i=l,6 
dummy 5 ( i ) =nc ( i ) 

450  continue 

if (debug. eq. . true. )  then 
write (15, *) 'guess 5' ,L 
call  flush(15) 
endif 

call  tenmatprod(Ce, strainc,sigec) 
call  rowcolumnprod(nc, sigec, dlam) 
dlam=dlam/L 
999  return 
end 

subroutine  yieldtest  (mhs,  sig,  sigyl,  f  ,yc,yf ) 
real*8  mhs(6,6) ,sig(6) ,dummy(6) , sigyl, f,yf 
logical  yc 
yc=. false. 

call  tenmatprod (MHS, sig, dummy) 
call  rowcolumnprod (sig, dummy, yf) 
yf =yf / ( (1 . 0-f ) *sigyl) 
yf =yf-sigyl 
if  (yf.gt.0.0)  then 
yc= . true . 
else 

yc= . false . 
endif 
return 
end 

This  routine  expresses  4th  order  tensors  as  6x6  matrices 
in  the  Voigt  notation 
subroutine  ten2matrix(ll, 1) 
realms  11 (3, 3, 3, 3) , 1 (6, 6) 

1(1, 1)=11(1, 1,1,1) 

1(1,2)=11{1,1,2,2) 

1{1,3)=11(1,1,3,3) 
l{l,4)=dsqrt(2.0d0) *11(1,1,1,2) 
l(l,5)=dsqrt(2.0d0) *11(1,1,2,3) 
l(l,6)=dsqrt(2.0d0) *11(1,1,3,1) 

1(2,1)=11(2,2,1,1) 

1(2,2)=11(2,2,2,2) 

1(2,3)=11(2,2,3,3) 
l(2,4)=dsqrt(2.0d0) *11(2,2,1,2) 

1(2,5) =dsqr t (2.0d0)*ll(2,2,2,3) 
l(2,6)=dsqrt(2.0d0) *11(2,2,3,1) 

1(3,1)=11(3,3,1,1) 

1(3,2)=11(3,3,2,2) 

1(3,3)=11(3,3,3,3) 
l(3,4)=dsqrt(2.0d0) *11(3,3,1,2) 
l(3,5)=dsqrt(2.0d0) *11(3,3,2,3) 
l(3,6)=dsqrt(2.0d0) *11(3,3,3,1) 
l(4,l)=dsqrt(2.0d0) *11(1,2,1,1) 
l(4,2)=dsqrt(2.0d0)*ll(l,2,2,2) 

1(4,3) =dsqr t (2.0d0)*ll(l,2,3,3) 

1(4,4)=(2.0)*11(1,2,1,2) 

1(4,5)=(2.0)*11(1,2,2,3) 

1(4,6)=(2.0)*11(1,2,3,1) 

1 (5 , 1) =dsqrt (2 , OdO) *11(2,3,1,1) 
l(5,2)=dsqrt(2.0d0) *11(2,3,2,2) 

1(5,3) =dsqr t (2.0d0)*ll(2,3,3,3) 

1(5,4)=(2.0)*11(2,3,1,2) 

1(5,5)=(2.0)*11(2,3,2,3) 

1(5,6)=(2.0)*11(2,3,3,1) 

1 (6,l)=dsqrt (2.0d0) *11(1,3,1,1) 
l(6,2)=dsqrt(2.0d0) *11(1,3,2,2) 
l(6,3)=dsqrt(2.0d0) *11(1,3,3,3) 

1(6,4)=(2.0)*11(1,3,1,2) 

1(6,5)=(2.0)*11(1,3,2,3) 


o  o 


1(6,6)=(2.0)*11(1,3,3,1) 

return 

end 


13 

12 

11 
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—  This  is  to  rotate  fourth  order  tensors 
subroutine  rot4order (11, q, Inew) 
real*8  11 (3.3 , 3, 3) ,q(3 , 3) , e (9) , f (27) , d(3) 
real*8  lnew(3,3,3,3) 
integer  i,j,k,l 
do  10  i=l,3 
do  11  j=l,3 
do  12  k=l,3 

f  (l)=q(l!l)*ll(l,l/l.l)+<3(l»2)*ll(l,l,l,2)+q(l,3)*ll(l,l,l,3) 
f (2)=q(l,l)*ll(l,l,2,l)+q(l,2)*ll(l,l,2,2)+q(l,3)*ll(l,l,2,3) 
f (3)=q(l,l)*ll(l,l,3,l)+q(l,2)*ll(l,l,3,2)+q(l,3)*ll(l,l,3,3) 
f (4)=q(l,l)*ll(l,2,l,l)+q(l,2)*ll(l,2,l,2)+q(l,3)*ll(l,2,l,3 
f ( 5) =q (If  1) *11 (1.2,2/ 1) +q (1 ,2)*ll(lf2,2f2) +q (1 , 3) *11 (1. 2,2, 3) 
f (6)=q(l,l)*ll(l,2,3,l)+q(l,2)*ll(l,2,3,2)+q(l,3)*ll(l,2,3,3) 
f (7) =q(l,  1) *11 (1,3 , 1,1) +q(l, 2) *11 (1,3, 1,2) +q(l, 3) *11 (1, 3, 1, 3) 
f (8)=q(l,l)*ll(l,3,2,l)+q(l,2)*ll(l,3,2,2)+q(l,3)*ll(l,3,2,3) 
f (9) =q (1,1) *11 (1,3, 3,1) +q( 1,2) *11 (1,3, 3,2) +q (1,3) *11 (1,3, 3, 3) 
f(10)=q(l,l)*ll(2,l,l,l)+q(l,2)*ll(2,l,l,2)+q(l,3)*ll(2,l,l,3) 
f (ll)=q(l,l)*ll(2,l,2,l)+q(l,2)*ll(2,l,2,2)+q(l,3)*ll(2,l,2,3) 
f (12)=q(l,l)*ll(2,l,3,l)+q(l,2)*ll(2,l,3,2)+q(l,3)*ll(2,l,3,3) 
f (13)=q(l,l)*ll(2,2,l,l)+q(l,2)*ll(2,2,l,2)+q(l,3)*ll(2,2,l,3) 

f ( 14 ) =q (1,1)*11(2,2,2,1) +q (1,2) *11 (2,2,2, 2) +q (1,3) *11(2, 2, 2, 3) 
f (15) =q(l, 1) *11 (2,2,3,l)+q(l,2)*ll(2,2,3,2)+q(l,3)*ll(2,2,3,3) 
f (16)=q(l,l)*ll(2,3,l,l)+q(l,2)*ll(2,3,l,2)+q(l,3)*ll(2,3,l,3) 
f (17) =q (1,1) *11 (2, 3,2,1 )+q( 1,2) *11 (2, 3,2,2 )+q( 1,3) *11 (2, 3, 2, 3) 
f (18)=q(l,l)*ll(2,3,3,l)+q(l,2)*ll(2,3.3,2)+q(l,3)*ll(2,3,3,3) 
f (19) =q(l, 1) *11 (3,l,l,l)+q(l,2)*ll(3,l,l,2)+q(l,3)*ll(3,l,l,3) 
f (20) =q (1,1) *11 (3, 1,2,1 )+q( 1,2) *11 (3, 1,2,2 )+q( 1,3) *11 (3, 1,2, 3) 
f (21)=q(l,l)*ll(3,l,3,l)+q(l,2)*ll(3,l,3,2)+q(l,3)*ll(3,l,3,3) 
f (22)=q(l,l)*ll(3,2,l,l)+q(l,2)*ll(3,2,l,2)+q(l,3)*ll(3,2,l,3) 
f (23)=q(l,l)*ll(3,2,2,l)+q(l,2)*ll(3,2,2,2)+q(l,3)*ll(3,2,2,3) 
f (24)=q(l,l)*ll(3,2,3,l)+q(l,2)*ll(3,2,3,2)+q(l,3)*ll(3,2,3,3) 
f (25)=q(l,l)*ll(3,3,l,l)+q(l,2)*ll(3,3,l,2)+q(l,3)*ll(3,3,l,3) 
f (26)=q(l,l)*ll(3,3,2,l)+q(l,2)*ll(3,3,2,2)+q(l,3)*ll(3,3,2,3) 
f (27)=q(l,l)*ll(3,3,3,l)+q(l,2)*ll(3,3,3,2)+q(l,3)*ll(3,3,3,3) 
e(l)=q(k,l)*f (l)+q(k,2)*f (2)+q(k,3)*f (3) 
e(2)=q(k,l)*f (4)+q(k,2)*f (5)+q(k,3)*f (6) 
e(3)=q(k,l)*f (7)+q(k,2)*f (8)+q(k,3)*f (9) 
e(4)=q(k,l)*f (10)+q(k,2)*f (ll)+q(k,3)*f (12) 
e(5)=q(k,l) *f (13) +q(k, 2) *f (14) +q(k, 3) *f (15) 
e(6)=q(k,l) *f (16) +q(k, 2) *f (17) +q(k, 3) *f (18) 
e(7) =q(k, 1) *f (19)+q(k,2) *f (20)+q(k,3) *f (21) 
e(8)=q(k,l)*f (22)+q(k,2)*f (23)+q(k,3)*f (24) 
e(9)=q(k,l) *f (25) +q(k, 2) *f (26) +q(k, 3) *f (27) 
d(l)=q(j,l) *e(l)+q(j,2) *e(2)+q(j,3) *e(3) 
d(2)=q(j,l)*e(4)+q(j,2)*e(5)+q(j,3)*e(6) 
d(3)=q(j,l)*e(7)+q(j,2) *e(8)+q(j,3)*e(9) 
lnew(i, j,k,l)=q(i,l)*d(l)+q(i,2)*d(2)+q(i,3)*d(3) 
continue 
continue 
continue 
continue 
return 
end 


This  routine  expresses  6x6  matrices  back  in  tensorial (4th  order)  form. 
The  matrix  is  in  Voigt  Notation, 
subroutine  mat2tensor (1, 11) 
real*8  1 (6, 6) , 11 (3, 3 , 3, 3 ) 

11(1,1,1,1)=1(1,1) 

11(1,1,2,2)=1(1,2) 

11(1,1,3,3)=1(1,3) 

ll(l,l,l,2)=l(l,4)/dsqrt(2.0d0) 

ll(l,l,2,3)=l(l,5)/dsqrt(2.0d0) 

ll(l,l,3,l)=l(l,6)/dsqrt(2.0d0) 

11(2,2,1,1)=1(2,1) 

11(2,2,2,2)=1(2,2) 

11(2,2,3,3)=1(2,3) 

ll(2,2,l,2)=l(2,4)/dsqrt(2.0d0) 

ll(2,2,2,3)=l(2,5)/dsqrt(2.0d0) 

ll(2,2,3,l)=l(2,6)/dsqrt(2.0d0) 

11(3,3,1,1)=1(3,1) 

11(3,3,2,2)=1(3,2) 

11(3,3,3,3)=1(3,3) 

ll(3,3,l,2)=l(3,4)/dsqrt(2.0d0) 

ll(3,3,2,3)=l(3,5)/dsqrt(2.0d0) 


11 (3 , 3 , 3 , 1) =1 (3 , 6) /dsqrt (2 . OdO) 

11 (1,2, 1,1) =1(4, 1) /dsqrt (2. OdO) 

11 (1,2, 2, 2) =1(4, 2) /dsqrt (2. OdO) 

11(1,2,3,3)=! (4, 3) /dsqrt (2. OdO) 

ll(l,2,l,2)=l(4,4)/2.0 

ll(l,2,2,3)=l(4,5)/2.0 

ll(l,2,3,l)=l(4,6)/2.0 

11 (2, 3, 1,1) =1(5,1) /dsqrt (2. OdO) 

11 (2, 3, 2, 2) =1(5, 2) /dsqrt (2. OdO) 

11 (2, 3, 3, 3) =1(5, 3) /dsqrt (2. OdO) 

ll(2,3,l,2)=l(5,4)/2.0 

ll(2,3,2,3)=l(5,5)/2.0 

ll(2,3,3,l)=l(5,6)/2.0 

11 (3, 1,1,1) =1(6,1) /dsqrt (2. OdO) 

11  (3, 1,2, 2) =1(6, 2) /dsqrt (2. OdO) 

11 (3, 1,3, 3) =1(6, 3) /dsqrt (2. OdO) 
ll(3,l,l,2)=l(6,4)/2.0 
11(3,1, 2, 3)=l(6,5)/2.0 
11(3,1, 3, l)=l(6,6)/2.0 
11(1,1,2,1)=11(1,1,1,2) 
11(1,1,3,2)=11(1,1,2,3) 

11(1, 1,1,3)=11(1, 1,3,1) 
11(2,2,2,1)=11(2,2,1,2) 
11(2,2,3,2)=11(2,2,2,3) 
11(2,2,1,3)=11(2,2,3,1) 
11(3,3,2,1)=11(3,3,1,2) 
11(3,3,3,2)=11(3,3,2,3) 
11(3,3,1,3)=11(3,3,3,1) 
11(2,1,1,1)=11(1,2,1,1) 
11(2,1,2,2)=11(1,2,2,2) 
11(2,1,3,3)=11(1,2,3,3) 
11(2,1,1,2)=11(1,2,1,2) 
11(2,1,2,3)=11(1,2,2,3) 
11(2,1,3,1)=11(1,2,3,1) 
11(2,1,2,1)=11(1,2,1,2) 
11(2,1,3,2)=11(1,2,2,3) 
11(2,1,1,3)=11(1,2,3,1) 
11(1,2,2,1)=11(1,2,1,2) 
11(1,2,3,2)=11(1,2,2,3) 
11(1,2,1,3)=11(1,2,3,1) 
11(3,2,1,1)=11(2,3,1,1) 
11(3,2,2,2)=11(2,3,2,2) 
11(3,2,3,3)=11(2,3,3,3) 
11(3,2,1,2)=11(2,3,1,2) 
11(3,2,2,3)=11(2,3,2,3) 
11(3,2,3,1)=11(2,3,3,1) 
11(3,2,2,1)=11(2,3,1,2) 
11(3,2,3,2)=11(2,3,2,3) 
11(3,2,1,3)=11(2,3,3,1) 
11(2,3,2,1)=11(2,3,1,2) 
11(2,3,3,2)=11(2,3,2,3) 
11(2,3,1,3)=11(2,3,3,1) 

11(1, 3, 1,1)=11(3, 1,1,1) 
11(1,3,2,2)=11(3,1,2,2) 
11(1,3,3,3)=11(3,1,3,3) 
11(1,3,1,2)=11(3,1,1,2) 
11(1,3,2,3)=11(3,1,2,3) 
11(1,3,3,1)=11(3,1,3,1) 
11(1,3,2,1)=11(3,1,1,2) 
11(1,3,3,2)=11(3,1,2,3) 

11(1, 3, 1,3)=11(3, 1,3,1) 
11(3,1,2,1)=11(3,1,1,2) 
11(3,1,3,2)=11(3,1,2,3) 

11(3, i,l,3)=ll(3, 1,3,1) 

return 

end 

SUBROUTINE  ZBRAK  { f  unc ,  XI ,  X2  ,  N,  XBl ,  XB2  ,  NB) 
real *8  xbl,xb2,xl,x2,x,dx, fp, fc, func 
integer  n,nb#i,nbb 
logical  yield,  check, neg 
logical  zb 

common  /mkdata2 /yield, check, neg 

common  /zbl/zb 

external  func 

NBB=NB 

NB=0 

X=Xl 

DX=(X2-X1)/N 

FP=func(X) 

if  (neg. eq. .true. )  then 
write (15, *) 'returning  in  ZBRAK' 
call  flushdS) 
return 


o  o 


i 


endif 

DO  11  1=1, N 
X=X+DX 
FC=func (X) 

if  ( zb . eq . . true . )  then 

write (15,  *)  'dlaiti'  ,X,  'f '  ,  fc 
endif 

if  (neg.eq. . true. )  then 
write  (15,  *j 'returning  in  ZBRAK__2  ' 
call  flush(15) 
return 
endif 

IF(FC*FP.LT.O.)  THEN 
NB=NB+1 
XB1=X-DX 
X32s=X 

write  (15,  *)  'info'  ,xbl,xb2,  fp,  fc 
call  flush(15) 

ENDIF 
FP=FC 

IF(NBB.EQ.NB)  then 

write (15, *) 'info2' ,xbl,xb2, func(xbl) , func(xb2) 
RETURN 
endif 
11  CONTINUE 
RETURN 
END 


